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1  INTRODUCTION 

The  optimal  control  model  (0<  'M  I  based  on  the  assumption  that  a  human 
operator  estimates  the  state  of  the  controlled  system  and  develops  a  control 
strategy  which  minimizes  a  performance  index.  The  pioneering  work  of 
Kleinman,  Baron,  and  Levison  [1].  of  Bolt,  Beranek.  and  Newman  (BBN) 
used  this  basic  assumption  to  set  up  an  optimal  control  problem  which 
closely  agrees  with  experimental  tracking  data.  The  resulting  controller 
consists  of  a  Kalman  Bucy  Filter  (KBF).  a  linear  predictor,  and  a  set  of 
Linear  Quadratic  Regulator  (LQR)  gains. 

Practical  understanding  of  the  OCM  develops  through  use,  which  in  turn 
requires  computer  implementation  Though  implementations  exist,  notably 
PIREP  [2],  they  are  not  readily  available.  The  purpose  of  this  report  is  to 
describe  an  implementation  of  the  OCM  using  Program  CC  [3],  This  will 
greatly  increase  the  availability  and  understanding  ol  the  OCM,  both  here 
at  STI  and  elsewhere. 

The  OCM  has  changed,  but  not  dramatically,  since  its  introduction  in 
the  late  1900's.  The  seminal  reference  [l],  see  also  (4,5],  has  the  following 
features:  a  performance  index  which  uses  a  weighted  sum  of  mean  square 
error  and  control  rate  energy,  full  attention  noise  ratios  for  observation  and 
motor  noises,  and  an  iterative  solution  method  to  achieve  the  desired  noise 
ratios.  None  of  this  has  changed. 

Additional  features  and  numerous  applications  have  appeared  in  the 
years  following  the  OCMs  introduction,  e.g.  [6]— [1 0).  The  additional  features 
included  in  Program  C’C's  implementation  are  visual  indifference  thresholds 
and  fractional  attention  parameters.  Notably  absent  is  the  use  of  pseudo- 
noise  to  induce  low  frequency  phase  droop,  and  t  he  optimization  of  fractional 
attention  for  multi-input  problems. 

The  current  interest  is  to  understand  and  predict  human  operator  be¬ 
havior  in  multi-axis  tasks  and  in  divided  attention  situations.  In  particular, 
the  objective  is  to  predict  pilot  behavior  and  ratings  for  multi-axis  tasks, 
and  to  compare  these  predictions  against  experimental  data.  Previous  work 
[1 1]-  [1-1]  has  suggested  that  the  value  of  the  optimal  control  performance  in¬ 
dex  correlates  with  pilot  opinion  ratings  (POR’s)  such  as  the  Cooper-Harper 
scale. 

The  motivation  for  implementing  the  OCM  on  Program  CC  is  to  further 
develop  the  ability  to  predict  POR's.  A  secondary  objective  is  to  simplify 
the  implementation  and  broaden  the  availability  of  the  OCM.  This  report 
describes  the  implementation  of  the  OCM.  Volumes  1  and  2  detail  some  of 
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the  applications. 

The  implementation  consists  of  several  Program  C C  macros  and  one  user 
defined  command.  The  macros  will  only  work  with  Version  4  of  Program 
CC.  The  macros  are  used  in  sequence  to  ( 1)  create  a  state  space  model  for  Yc 
and  Yw ,  the  controlled  element  and  driving  noise  filter,  (2)  iterate  the  LQR 
problem  in  order  to  set  the  neuro-muscular  mode  or  modes,  (3)  set  up  the 
KBF/linear  predictor  problems.  (4)  iterate  the  KBF  problem  until  desired 
noise  ratios  are  achieved,  and  (5)  analyze  the  resulting  Yp  pilot  model.  The 
user  defined  command  is  called  from  within  the  macros  as  part  of  the  LQR 
and  KBF  iterations.  In  addition  to  the  macros,  all  of  the  existing  and 
powerful  capabilities  in  Program  CC  can  be  used  for  analyzing  the  resulting 
system  models  and  pilot  models. 

Section  2  contains  a  technical  description  of  the  OCM.  It  is  presented 
in  its  general  multivariable  form,  but  only  the  single-input  version  is  imple¬ 
mented.  The  description  brings  together  material  scattered  among  several 
references  [l,2],[4]-[9j,  and  adds  a  new  twist  to  the  treatment  of  delay  which 
results  in  state  space  and  transfer  function  models  for  the  pilot.  Section  3 
is  the  main  reference  for  the  Program  C'C  implementation;  providing  usage 
notes,  and  listing  in  Tables  1-3  the  macro  names,  macro  parameters,  and 
data  storage  locations.  Several  examples  are  presented  in  Section  4,  includ¬ 
ing  the  use  of  the  optimal  cost  to  predict  pilot  opinion  ratings.  Background 
information  on  Program  C’C  is  provided  in  Appendix  A,  and  the  macros  are 
listed  in  Appendix  B. 

Fortunately  it  is  not  necessary  to  completely  understand  the  technical 
details  of  the  OCM  in  order  to  use  it.  It  is  sufficient  to  use  the  OCM  as 
a  means  to  an  end:  a  long  handled  crank  which  results  in  a  )'p.  Assuming 
some  prior  knowledge  about  the  OCM  and  about  Program  CC.  the  minimum 
amount  of  addition  information  needed  to  operate  the  macros  is  contained 
in  Tables  1  and  2. 
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2  TECHNICAL  PROBLEM  STATEMENT 


A  complete  technical  description  of  the  optimal  control  model  is  presented 
in  this  section.  The  description  is  valid  for  the  multivariable  case,  though 
only  the  single-input  single-output  case  is  implemented.  The  equations  are 
summarized  in  the  block  diagram  of  Figure  1,  and  the  computational  flow 
is  summarized  in  Figure  2. 


2.1  Controlled  System 

The  controlled  system  is  modeled  using  state  space  equations.  The  driving 
noise  w(t)  will  typically  contain  dynamics,  which  are  included  together  with 
the  system  dynamics  in  the  A  matrix.  The  only  difference  between  the 
single-input  and  multi-input  versions  is  the  dimension  of  the  input  vectors. 


i(t )  =  Ax(  t)  -f  Bu{t)  +  Ew(t) 

2/(0  =  Cx(t)+Du(t) 
yP(t )  =  y(t  -  r)  +  vy(t  -  r) 

The  dimensions  and  definitions  of  the  vectors  are  as  follows.  The  notation 
Rn  indicates  a  vector  of  n  r-jal  numbers. 

r(t)  €  Rn  state  yp(t)  6  Rr  observation 

u(t)  G  Rm  input  </'(/)  €  Rm'  driving  noise 

y{t)  £  Rr  output  !’,(/)£  If  observation  noise 

The  noise  intensities  arc: 

E{  «•(  l)w(  /  —  rr )}  =  r„«f»(rr) 

E{ry{t)vv(l  -  a)}  = 


It  is  assumed  that  the  human  uses  both  errors  and  error  rates,  and 
therefore  the  output  y(t)  contains  both.  It  is  not  important  how  they  the 
outputs  are  arranged,  so  for  notational  convenience  the  errors  are  grouped 
together  and  listed  first: 


i/(0  = 


(!hd) 

Uoo 


Tim  transfer  function  Vf(s)  for  the  controlled  system  is  defined  to  be  from  the 
input  u  to  the  error  ye,  because  the  single- input  case  this  is  the  traditional 
definition.  The  Program  C’C  macros  described  in  Section  3  use  variations  of 
Vc(.s)  with  different  inputs  and  outputs. 
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Figure  1:  The  Optimal  Control  Model 


(i 


The  human  visual  <  1*  I  a  v  nl  -  >i'nniih  is  included  in  the  definition  of  the 
observation  yp(t).  TypnalU  lot  the  ()(  M  r  -  .2  is  use<i.  Note  that  this  t 
is  not  the  same  value  as  the  effective  delay,  r, .  which  is  based  on  the  unit 
magnitude  crossover  of  V,,V  .  Fxperimentally  determined  values  for  r,  range 
from  .15  to  .25  seconds  and  depends  on  the  amount  of  pilot  lead. 

The  driving  noise  intensity  l„  is  determined  as  part  of  the  experimen¬ 
tal  setup.  The  observation  noise  intensity  \y  is  automatically  determined 
according  to  an  implicit  equation  explained  in  Section  2.5. 

2.2  Performance  Index 

’flie  objective  of  the  OCM  is  to  determine  a  control  law  u(l)  =  f(yp(t)) 
which  minimizes  the  performance  index: 

.1  =  lim  E  j  —  /  x'Qs  -t-  u(lu)<lt 

I  —  s.  y  1  Jt) 

In  most  if  not  all  of  the  cases  in  w  hi  eh  t  he  ()( 'M  has  been  used,  the  weighting 
matrices  take  the  diagonal  form: 

Q  --  '"W 

Qy  -  diag^r/i . </•■} 

(>  -  d  i  a  g  { <i  i . <;  . } 

The  (/,' s  corresponding  to  the  errors  are  scaled  to  give  approximate'  eepial 
weight  to  standard  deviations,  and  the  e/,'s  corresponding  to  error  rates  are 
set  t  o  zero: 

1  /<  . r/2 

0  i  -  r/ 2+1 . r 

The  i/,'s  are  chosen  to  in  order  to  plaee  the  neuromuscular  dynamics,  as 
explained  m  the  next  subsection. 

1  he  ()CM  differs  from  tlm  standard  l  inear  Quadratic  Gaussian  (I.QG) 
opt  imal  tout  rol  problem  in  t  hree  respects:  (  I  |  t  lie  weighting  on  control  rate, 
i  2  the  inclusion  of  motor  noise  before  the  neuro muscular  dynamics,  and 
( .1  the  observation  delay  of  r  sooimK. 

2.3  LQR  Solution 

i  he  stochastic  Linear  Quadratic  Regulator  problem  is  solved  to  determine 
t  lie  full  st  at e  <  out  rol  weights  The  system  is  augmented  w  ith  a  preceding  in¬ 
tegrator,  which  results  in  an  equivalent  performance  index  without  a  control 


s 


rate. 


io{t) 

=  -d0J',)(  t  I 

B,jii  -  A(i  ic 

\ 

J 

=  Inn  E 

T—o c 

t  T7  /  (‘r'j^°Xo  +  tl’Gtl)  dtj 

where: 

-do  = 

(A  B V 
VO  o ) 

*-(!) 

*-(f) 

Qo  = 

(Q  o] 

V  o  o ) 

n  =  ii 

The  solution  is: 

A 

-  -  A  Jo 

T 

tj 

=  CrlB'0K0 

where  A‘0  is  the  unique  positive  semi-definite  solution  of  the  algebraic  Riccati 
equation  (ARE): 

0  =  Aq/vo  +  AqAo  +  Qo  —  A  o  lh)(<  1  BqBo 

The  solution  method  of  choice  for  the  ARE  is  to  ( 1 )  set  up  the  Hamilto¬ 
nian  matrix,  (2)  compute  its  Schur  decomposition.  (3)  order  the  eigenvalues 
of  'he  Hamiltonian  so  the  stable  eigenvalues  are  in  the  upper  left  hand 
block  7'n  of  upper  triangular  Schur  matrix  T,  and  then  (1)  compute  A’o 
using  the  partitioned  Schur  vectors  (i  which  span  the  stable  subspace.  This 
is  a  modified  Potter's  method,  due  to  haul).  Schur  vectors  are  computed  as 
an  intermediate  step  towards  eigenvectors,  the  latter  of  which  are  in  general 
numerically  ill-conditioned. 

(  -do  _  fl  n  IvAf'lw  f  1 2  A 

\  -Qt)  -do  /  V  ^  21  f  22  /  \  d  T-ll)  \^21  l  22  ) 

Bq  -  i  21^  1 1 ‘ 

The  control  gains  can  be  partitioned  to  separately  apply  to  the  state  r 
and  the  input  u: 

A  -  (  At  £2)(') 
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By  convent  ion  t  In'  feedback  weinht  on  u  is  closed  prior  to  t  lie  ot  her  feedbacks 
resulting  in  the  neuro-nniscular  dynamics.  Input  to  the  neuro-muscular 
dynamics  then  consists  of  the  remaining  feedbacks  plus  injected  motor  noise. 

T\ii  +  u  -  ua  +  rUo 

ua  =  -l.’r  o 


where: 


Tv  -  L~2 1 

L'  -  (  T71  L\  0) 

and  where*  the  motor  noise  intensity  is: 

E{rUo(/)i’Uo(/  -  er)}  =  Vu.,Ha) 

The  value  of  V’Uo  is  automatically  determined  as  described  in  Section  2.5. 

The  input  rate  weight  G  is  adjusted  so  as  to  achieve  neuro-  muscular 
dynamics  of  10  rad/scc.  In  the  single-input  case  a  simple  binary  search 
yields  the  correct  value  of  G.  In  the  case  of  decoupled  inputs  separate 
binary  searches  ran  be  used  for  each  input.  For  cases  more  complicated  no 
help  is  available. 


2.4  KDF  and  Linear  Predictor  Solutions 

The  states  are  not  available  for  feedback,  only  the  delayed  and  noise  cor¬ 
rupted  observation  yp(t).  The  optimal  solution  proceeds  by  using  a  KBF  to 
estimate  the  state  at  time  t  -  r.  and  then  using  a  linear  predictor  to  esti¬ 
mate  the  state  at  time  t.  Due  to  the  convention  used  for  the  neuro-inuscular 
dynamics,  the  following  augmented  system  is  used  for  the  KBF: 

Fi(0  =  .4|.ri(f)  +  B\lla  +  (ri 

yp(  t]  ~  C  lr  1  (  t  —  ~  )  +  ry(  t  -  T  ) 

where: 
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The  KHF  computes  p{t)  the  linear  mean  square  estimate  of  X\{t  -  r) 
based  on  observations  yp(rr)  for  n  <  /: 

p{t)  =  Aip(t)  +  // 1  [i/p(n  -  Cip(t)}  +  Bxua(t) 

where  the  filter  gain  is: 

//,  =  zlc[v-' 

and  where  Ej  is  the  unique  positive  semi-definite  solution  of  the  ARE: 

o  =  a,e,  +E1.1;  +  iv,  -E.rjr-'c.s, 

The  linear  predictor  updates  p(  t )  to  obtain  x  i  ( / ),  the  linear  mean-square 
estimate  of  X] (t)  based  on  observations  yp{n)  for  a  <  t.  Note  that  p(t-r)  ^ 
Xl(0- 

£(0  =  Alt(t)+  BXU'AD 

x,  =  S(t)  +  <y'hT[p(t)-Z(t  -T)} 

The  LQR  control  weights  are  applied  to  i\(t),  again  used  according  to 
the  convention  for  the  neuro-muscular  dynamics: 

u»(t)  =  -rh(t) 


2.5  Noise  Ratios,  Indifference  Thresholds,  and  FVactional 
Attentions 

The  observation  noise  intensity  lv  ami  the  motor  noise  intensity  l’Uu  are 
assumed  to  Ire  diagonal  matrices: 


Vy  =  dia^j  ,VI . 

VUa  =  diag{  rllai . . . . 
where  each  of  the  diagonal  elements  satisfy: 

y'  f,  E{try,.r,)* 


*  r««,r.  ) 


i  =  1 . r 


v 


Ua, 


P^a, 


na„ 


i  =  1 . in 


(1) 

(2) 


1 1 


and  where  the  individual  terms  in  { 1 .2  ■  are  defined: 

Pyx'  Pua,  ~  noiso  ratios 

<r y, .  aUa  —  standard  deviations  of  y,,  uai 
f,  =  frartion  attentions 
T,  —  indifference  thresholds 

E<c„,T,)  = 

Equations  (1.2)  for  the  no'se  intensities  vy>  "s  and  rUo  ’s  are  implicit, 
because  the  standard  deviations  tfVt"s  and  oUa  's  depend  on  the  noise  inten¬ 
sities.  An  iteration  is  required  in  order  to  compute  the  noise  intensities. 
The  iteration  takes  t he  following  form:  (1)  guess  initial  values  for  the  noise 
intensities.  (2)  compute  the  KBF  gain  followed  by  the  standard  devi¬ 
ations,  (3)  if  equations  (1.2)  are  satisfied  within  a  specified  tolerance  then 
stop,  else  use  equations  (1.2)  to  compute  new  values  for  the  vVl  "s  and  vUa  's, 
and  then  repeat. 

The  observation  noise  ratios  pVi  have  been  experimentally  determined 
to  be  .01.  Usually  this  ratio  is  expressed  using  power  dB  101og10,  or  in  this 
case  -20  dB. 

The  motor  noise  ratios  can  be  measured  in  principle,  but  there  has  not 
been  any  found  in  practice  (when  dither,  stick  pumping,  etc  is  excluded). 
Various  hypotheses  have  been  proposed  for  the  origin  of  motor  noise,  but  it 
is  more  likely  that  motor  noise  is  a  mathematical  fiction  used  to  achieve  a 
solution.  Numbers  used  for  the  />Uo  's  vary  from  .001  to  .05,  some  researchers 
preferring  the  smaller  values  and  some  the  larger.  A  “safe"  value  seems  to 
be  .01. 

The  output  ij,  must  move  past  an  indifference  threshold  T,  before  a 
change  is  perceived.  This  threshold  is  modeled  by  a  Gaussian  describing 
function.  The  input  of  the  threshold  at  a  particular  time  is  normally  dis¬ 
tributed  wiui  mean  zero  and  standard  deviation  oyi.  and  the  output  is  nor¬ 
mally  distributed  with  mean  zero  and  standard  deviation  ayJ E{oyi.T,), 
where: 

E{oy>,T,)  =  erfc  ^ =  >  T,) 

where  x  is  a  sample  from  a  normal  distribution  with  mean  zero  and  standard 
deviation  <ryi. 
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When  the  indifference  threshold  is  7,-0  then  £=1.  and  it  is  always  the 
case  that  0<£'<1.  Visual  perception  thresholds  have  experimentally  been 
determined  to  be  .05  deg  and  .1  deg/sec.  The  value  of  7,  used  in  the  OCM 
depends  on  how  far  the  display  is  from  the  human. 

The  observation  noise  is  assumed  to  be  inversely  proportion  to  the  frac¬ 
tional  attention  /,,  e.g.  the  observation  noise  intensity  vVi  doubles  if  only 
50%  attention  is  being  paid  to  output  y,.  It  does  not  mat  ter  how  the  remain¬ 
ing  attention  is  allocated,  whether  for  other  control  tasks  or  for  non-control 
tasks  such  as  communication.  The  usual  assumption  is  that  the  same  atten¬ 
tion  is  attributed  to  the  error  and  error  rate  in  a  given  axis,  therefore: 

£/■  =  > 


Or  in  terms  of  the  outputs: 

E/-2 

i=i 

Some  researchers  allow  /,  to  differ  for  error  and  error  rate,  and/or  require 
the  total  across  outputs  to  be  1.  In  the  single-input  (single-axis)  case  simply 
use  /i= /2  =  1.  In  the  multi-axis  case,  if  the  tasks  are  of  comparable  difficulty 
then  apply  equal  attention  to  each  axis,  otherwise  optimize  over  the  /,  so  as 
to  minimize  the  performance  index  ./.  Program  C'C  allows  arbitrary  input 
of  /,  and  docs  not  check  the  sum,  but  does  not  optimize  over  /,. 

Pseudo  motor  noise  is  a  scheme  where  large  values  of  pUa  are  used  to 
compute  the  KBF  gains,  and  then  smaller  values  are  used  to  compute  the 
standard  deviations  oy> \s  and  oUa  's.  Different  equations  then  those  pre¬ 
sented  in  Section  2.5  are  needed  in  order  the  compute  the  standard  devia¬ 
tions,  due  to  the  fact  that  the  state  estimation  errors  are  no  longer  uncor¬ 
related  with  the  state  estimates.  The  use  of  pseudo  motor  noise  has  been 
proposed  to  account  for  phase  droop  and  motion  cues.  Pseudo  motor  noise 
is  not.  however,  included  as  part  of  the  Program  CC  implementation,  in 
other  words  the  same  values  for  pUa<  must  be  used  for  all  calculations.  This 
is  no  big  loss.  Pseudo  motor  noise  is  a  historical  oddity  for  which  no  general 
guidelines  were  developed  and  which  never  lived  up  to  its  expectations. 

The  7r  in  equations  (1,2)  was  the  source  of  confusion  to  the  author.  Even 
though  others  are  not  likely  to  be  similarly  fooled.  Til  take  the  liberty  to  ex¬ 
plain.  Neglect  for  this  argument  the  attentional  fraction  and  the  indifference 
threshold,  leaving: 

Vy  =  Py^ay 
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The  mean  square  value  of  ih«-  out  pm  is  defined  using  t he  autocorrelation 
function: 

i  -  «y(0) 

where  the  autocorrelation  is  defined: 

HyiT)  =  E \y(t)y(t  -  t)} 

and  where  the  autocorrelation  is  related  to  the  power  spectral  density  using 
the  Fourier  transform  identities: 

Sy(ju)  =  \  Ry(T)e~JUJTdT 

7—  x 

Ry(r)  =  Sy(juj)cJ^Tdu; 

The  noise  ratio  py.  in  terms  of  the  power  spectral  density  is: 

Vy 

P'J  SyU~')<L: 

If.  however,  the  integration  of  Sy(j\j)  is  defined  only  over  the  positive  fre¬ 
quencies  then  the  noise  ratio  is: 

=  v> 

Py  fo^  Sy(j<j)(U 

hence  py  can  be  defined  as  ratio  of  the  input  noise  intensity  to  the  positive 
frequency  power  density. 

The  OCM  literature  over  the  years  defines  the  noise  ratio  py  using  posi¬ 
tive  frequencies  (where  no  tr  is  needed),  but  uses  calculations  in  terms  of 
(where  t  is  needed).  To  make  a  long  story  short,  use  n  in  the  calculations. 

2.6  Performance  Measures 

The  time  domain  performance  measures  are  the  mean  square  errors  and  the 
optimal  cost.  The  mean  square  for  the  augmented  state  is: 

A  =  E{*,x/,}  =  e>,'rEie'4,«r+  /  e^ir.e^dcr 

J o 
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Define  the  estimation  error  >  \  —  X\  s^  I  lien  X  can  be  decomposed  into 
X  =  E  +  A  ,  where: 


E{ej6i}  =  e','rAie‘,i  - 

E{i,i;}  =  r~{A'  HJ 

Jl\ 


/  e4'/TirIe-4>ffd<T 
”  *>A'T  H  xYylI\eA\T  q{A'~b'l">'°  do 


The  augmented  state  is  j'J  —  (  s'  it' ;.  hence  the  mean  square  values  for  x 
and  u  are  respectively  the  A’u  and  V22  blocks.  The  mean  square  output 
and  neuro-muscular  input  are: 


V  =  E  {yy}  =  <:\XC[ 

ua  -  =  r.vi 


The  mean  square  of  u  is  infinite  because  the  white  motor  noise  feeds  directly 
into  it.  Using  instead  the  estimated  value  of  it: 

Udot  =  E{  it  u}  =  L  X  £'  +  (  0  L2 )  E  (  0  Li)' 


Finally,  the  optimal  cost  is: 

-/  =  lim  E  l  rr  [  ( !}  Q  v  U  +  it'd  it)  <lt  \ 
r— ^  /  Jo  j 

=  Trace  [E {y<j'}Qy}  +  Trace  [E{«i/}6'] 

=  Trace[VQv]  +  TracefC,^//'] 


which  in  the  single-input  case  is: 

=  (ri  + 

The  types  of  calculations  needed  for  these  performance  measures  are  ma¬ 
trix  exponentials  and  Lyapunov  equations.  The  preferred  method  for  matrix 
exponentials  is  the  scaled  Fade  approximation,  (with  the  scaled  Taylor  se¬ 
ries  a  close  second  and  the  unsealed  Taylor  series  absolutely  not  to  be  used). 
The  scaling  is: 

A  =  A\t /a,  where  o  =  22 

and  where  m  is  the  smallest  integer  >0  such  that  n  is  greater  than  the 
square  root  of  the  sum  of  squares  of  the  elements  of  Air.  This  scaling 


15 


approximately  keeps  the  eigenvalues  of  Ax  r  less  than  1.  If  F  =  exp(/l)  and 
F  =  exp(AiT),  then  rescaling  is  accomplished  simply  by  F  —  F2  ,  which 
requires  only  m  matrix  multiplications. 

A  block  matrix  is  used  to  compute  the  integral  of  a  matrix  exponential: 


exp 


where: 


F,  =  e"4‘r 
F2  =  er>T 


Hence: 


G'i  =  I*  o-A'lT-°)WlcA>ed(T  =  e--4ir 

Jo  Jo 

E  =  e'4,rE,ei4!T  +  (*  oA'0\\\oA'ado  =  F2%F2  +  F2'G, 
Jo 


The  Lyapunov  equation  is: 


0  =  PQ  +  QP  +  R 


where  R  is  symmetric.  The  solution  Q  exists  and  is  uni<|iie  if  the  eigenvalues 
of  P  and  R  do  not  coincide.  This  condition  is  automatically  satisfied  if  P  is 
stable  and  R> 0.  If  P  is  stable  then  the  Lyapunov  solution  Q  satisfies  the 
integral  identity: 

Q=  r  c/v  Rcp'*do 
Jo 

Hence  the  mean  square  value  A  is  computed  using  the  Lyapunov  equation 
with: 


P  =  ,»,  -  I  UP 
It  =  e  *1  T II 1 1  V//Je ■'> r 


The  preferred  solution  method  is  Bartel-Stewart,  using  the  Schur  decompo¬ 
sition  for  P. 


2.7  Frequency  Responses 

The  frequency  responses  of  interest  are  >'c(.s).  Yp(s ),  and  the  remnant 
$nn((6).  Here  it  is  explained  how  to  compute  the  frequency  response  data, 
and  a  state  space/transfer  function  approximation  for  }'p(s). 

The  pilot  model  Yp(s)  is  composed  of  a  visual  delay,  KBF,  linear  predic¬ 
tor,  and  neuro-muscular  mode. 

P  =  (.4i  -  H\C\)p  -t  Hxy(t  -  t)  -f  Bxua(t  -  r)  +  H\i’y(t  -  r) 

i  -  A\£  +  Biua 

ii  =  ^(/)  +  e-4ir[p(t)-^(/ -  r)\ 

ua  =  —  L'x  i 

u  =  -L2U  +  L?ua  +  L-ji'u 

It  is  more  convenient  to  replace  t he  differential  equation  for  £  with  another 
for  ix: 


i  1  =  -eA'TIIxCxp  +  (Ax  -  BxLm)x  1  +  eA>T Hiy  +  eAtT  Hx vy 


Move  the  visual  delay  from  y  over  to  «a.  which  is  equivalent  from  an  input/ 
output  point  of  view.  Combine  p  and  ii  into  an  augmented  vector  x 2. 
Change  the  input  matrices  for  y  to  obtain  an  equivalent  system  driven  just 
by  ye.  The  result  is: 

■i'2  =  A 2X2  +  fh'Jt  +  £2 Unit  -  T)  +  B2 Vy 

H.i  =  ('2*2  +  Ihilr 

a  =  -L2U  +  l<2 Mu{/  -  r)  +  l,2fu 


where: 

f  A\  -  U\C\  0  \ 

2  \-eA'rHxCx  AX-BXL’J 


lh=  (eA/J  =((52)l  (/?2):) 


x2  = 


C’2  =  (0  -L*  ) 


B2  -  (#2)1  +  ^2  -  Of  #2)2 


The  creation  of  fi2  and  Z)2  is  based  on  the  identity: 


s(sl  —  /42)  1  —  /  +  (s/  —  .4 2 )  */l2 
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Frequency  data  points  can  be  exactly  calculated  by  using  l In*  Laplace 
transform  e-sr  of  tin*  delay.  I  lie  following  steps  result  in  calculation  for 
ip(s)  and  ^n„e(s).  Define  the  following  transfer  functions,  located  as  shown 
in  Figure  3. 


H(s) 

H3(s) 

F(s) 

L  mi  $) 

Vc(.s) 

V'u..(«) 


(  Hl(fi)  ffi(fi)  )  =  f-iisl  -  /I2)  ’#2 
H,(s)  +  sll2(s)  =  C2(sl  -  A2)-'B,  +  b2 
C2(sl  -  a2)~x  e2 

(si-  L2r'L2 
(C)i(sl  -  A)~'B 
(Ch(sI-A)-'E 


where  C  - 


Tne  transfer  function  from  ye(s)  to  »„(*)  is  the  solution  of  the  implicit 
equation: 

«a(s)  =  //3(  •“)!/<(•$)  +  E(s)v~,Tua(s) 

which  is: 

ua(s 1  =  [/-  /r(s)c-*r]~l  Il:i(s)yf(s) 


The  pilot  model  Vp(s)  is  the  transfer  function  from  y((s)  to  u(s): 


Yp(s)  =  Lm(s)Gf(s)Ih(s) 

where  6’/(,s)  =  e_,r  [/  -  F(s)e-5r]  1 


In  order  to  determine  the  remnant,  use  Figure  3  to  determine  the  transfer 
functions  from  the  noise  inputs  to  the  error: 


ih (*)  =  (/  -  YpYc)  '  (  YeLmG/ll 

' - S/ - 


The  reinant  $„„,(«)  is  the  spectral  density  resulting  from  the  i\(s)  and 
i'ua(s)  noise  sources: 


=  Gci(s)d\ag{Vy.  VUa.  Q)G'c,(-s) 

Up  to  this  point  the  visual  delay  of  r  seconds  has  been  exactly  modeled. 
It  is  not,  however,  possible  to  obtain  a  state  space  description  of  Yp  unless 
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an  approximation  is  used  (o'-  'n<-  Obtaining  a  state  i space  realization 

is  quite  informative,  because  b>  converting  to  transfer  function  for  }’p(s), 
the  poles  and  zeros  can  be  examined.  \  state  space  realization  of  the  Fade 
approximation  of  a  delay  is: 

Xp(t)  -  .-IpjTplM  +•  HpUri(t) 

up(t)  -  CBjr„{t)  -f  Dpua(t) 

whore  up{t)  s;  ua(t  -  r) 

For  example,  a  2nd  order  Fade  approximation  for  the  single-input-single¬ 
output  case  is: 


Up{t)  =  (0  -12/7- )*„(/)  +  (  1  )Ua(t) 


The  2nd  order  case  usually  suffices,  though  others  may  of  course  be  tried. 
Substitute  the  Pade  approximation  and  group  all  of  the  differential  equation 
together  to  obtain  the  following.  The  input  ij  can  be  compressed  as  described 
earlier  in  order  to  obtain  an  equivalent  stale  space  system  driven  just  by  ye. 


Frequency  data  for  V’p(s)  can  be  obtained  directly  from  the  above  state 
space  system.  In  order  to  obtain  the  remnant  a  state  space  description  of 
the  closed  loop  system  must  be  obtained,  but  rather  than  attempt  to  write 
it  out  in  terms  of  all  of  the  elementary  matrices,  it  is  best  to  let  a  program 
such  as  ('('  do  the  required  algebra.  The  closed  loop  system  Gc[(s)  which 
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will  be  of  the  form: 


j-4 
Ve 

As  before,  the  remnant  is: 

^nnt(s)  =  6',/(*)<liag{lV  ^  tiu'  0}G'cl(-s) 

Several  quantities  of  interest  can  be  obtained  from  Bode  plots  of  V'p(.s) 
and  Yp(s)Yc(s).  including: 

w\-  =  unit  magnitude  crossover  frequency 

Qp.\i  =  phase  margin 

Te  =  (~/2  -  Opm  )/wV  =  effective  delay 

Op  =  -ArgppOw'c)]  +  w.r,  +  arctan[w-crv] 
pilot  phase  compensation 

These  are  best  obtained  by  using  a  cursor  to  read  magnitudes  and  phases 
from  Bode  plots. 
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3  OCM  IMPLEMENTATION 

3.1  Overview 


The  human  optimal  control  model  is  implemented  as  a  series  of  macros, 
two  of  which  call  a  user  defined  command1  The  following  functions  are  per¬ 
formed: 

1.  setup  a  state  space  model  of  the  controlled  system 

2.  solve  the  LQR  problem 

3.  setup  the  KBF  problem 

4.  compute  KBF  iterations  and  performance  measures 

5.  compute  state  space,  transfer  function,  and/or  frequency  response 
models  of  the  human  operator 

Both  the  LQR  and  KBF  problems  involve  iterations,  which  are  auto¬ 
matically  controlled  by  the  user  defined  command.  After  convergence  is 
obtained,  all  of  the  regular  features  of  Program  CC  are  then  available  for 
analysis  of  the  results.  The  types  of  analysis  include  frequency  responses 
and  simulations.  With  the  ability  to  obtain  transfer  function  models  of  the 
human  operator,  the  types  of  analysis  can  be  extended  to  include  the  root 
locus,  and  more  importantly  to  include  low  order  transfer  function  approxi¬ 
mations.  [t  is  these  low  order  approximations  which  clarify  the  relationship 
between  the  optimal  control  and  classical  pilot  models. 

I  lie  information  for  the  macros  is  summarized  in  Tables  1  through  3. 
Only  the  SISO  case  is  implemented. 

1.  Table  1:  OCM  Macros 

2.  Table  2:  Input  Parameters 

3.  Table  3:  Output  Parameters 


*A  niacrn  is  an  indirect  command  file,  whereas  a  u ft  r  dtjintd  command  is  a  compiled 
BASIC  program  which  can  he  chained  to  Program  CC 
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Table  1:  OCM  Macros 


_ I 

1.  OCMYCl.G,.C,,(Vc)„„ 

Create  state  space  (}'  )„//  for  output  driving  noise. 
ye  =  6’,(s)u  +  Gj(s)ir 

2.  OCMYC2.  6',,  Gj,  {Yc)an 

Create  state  s/>ace  ( )’c)ait  for  input  driving  noise 

Ue  =  G,{*)[u  +  ( .s )  u-] 

3.  OCMYC3,  G’„6V  6V  (>;-)„» 

Create  state  s/xtee  for  gene  ml  driving  noise 

y?  =  G,(s)[Gj(s)u  + 

■1.  OCMALL.  (V'c)a„ 

Calls  OCMLQR,  OCMSFTBB.  and  OCMKBF 

5.  OC  MLQR.  0 c)aih  (h  9'  Tv*  thrrsk 

Automatic  iteration  of  g  to  achieve  desired  t\ 

().  OCM  St  TIP.  ru„  r,  Pyx.  py2.  pUa.  lVl.  lVJ.  V„a.  r#1,  Tj,;.  / 

Setup  KBF  and  linear  predictor  problem 

I 

|  7.  OCM  KB  Y.dBlhr,sk 

j  Automatic  iteration  of  V  s  /or  KBF  problem 

I 

v.  OCMPILOT,  r.  Parle  order 

Create  state  si>ace  models  ofYp(s)  and  }'.;(*) 

!  ».  OCMG.  G’,.  6’ 

i 

Create  tfs  ) ;)  rmr/  I  r  fmm  state  space  mrxlels 

in.  OCMFREQl,^.,..^,,*.  #pts 

Blot  of  Yp,  Yr.  YPYC,  and  using  state  space  results 

j  11.  OOMFRF.Q2.  r .  *jUj,v .  #pts 

!  Blot  of  Yr,  Yc.  YpYc.  and$nrlr  using  e-sr 


Table  2:  Input  Parameters 


.Vote:  ( Yc)aii  ,s  n  stuff  space  quadruple.  Do  not  use  names 
which  overlap  those  used  by  the  macros:  P  — *  P30  and 
P$ oo  — ■ *  D 504 -  T 40  has  become  a  standard  name  to  use 
for  (Yc)all-  The  parameters  C,,  Gj.  and  G k  are  transfer 
functions,  for  which  any  names  can  be  used.  All  of  the 
remaining  parameters  are  entered  as  numbers. 

=  controlled  system  and  driving  noise  dynamics,  where: 
input  =  ^  j  output  - 

<?,  <j  -  quadratic  weights,  where: 

fo  (< iyl  +  dt 

t\  -  desired  neuro  muscular  time  constant 
Tthnsh  -  desired  accuracy 
l'„.  =  driving  noise  intensity 
r  =  visual  delay 

Pm  -  I'm-  P<‘ „  =  noise  ratios  (entered  in  dli's) 

\  ,,, .  l  Vi,  r„a  =  observation  and  motor  noise  intensities 
f  \ .  7j  =  indifference  thresholds  (0  =  no  threshold) 

/  —  fractional  at  ton' ion  (1  —  full  attention) 
dDthresh  =  desired  dB  accuracy  of  noise  ratios 
Pade  order  =  order  of  e_ir  Pade  approximation 
.  *nigh.  #pts  =  Bode  plot  axis  limits 


J  =  lim  e!  4 


T-cc 


T  . 
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Table  3:  Output  Parameters,  Single-Input  Systems 


Created  by  OCMLQR: 

I’  =  (Yc)all  Pj,  =  g 

P\  =  not  used  P.\  —  L  =  (  l.\  Li) 

Pi  =  <1 

(noted  by  OCMSETVP: 

>  P\o  =  FA  \  E1 

=  (  //j’Z-i  0  )  T|i  =  A’Z?/-  iteration  counter 

,n  =  (si  -  L2)~l  L 2  P\2  =  (Py>  Py2  />u,  7y,  7y,  /  r 

1  (  -^  /  —  T  1  )  “  1  /i  1  /'n  =  (  V«1  '  «  '  II  a  ) 


Created  by  OCMKIil : 

Pj  9  =  V  =  /;  +  .V 
11 }  ,  7>20  =  (  <  <  <,  ) 

lt>  e',lf,U’i o^Cela  P2]  =  not  used 

‘A>  rT1e/,i r  (£=/»„;+/», 7 )  P22=  not  used 
V  =  (  '/"y,  7  ) 

('noted  by  OCM PILOT  (state  s /xi re  models): 

Pn  =  vp 

/b.s  =  V  /  =  C.|( .<•/  -  .1 ,  )“ 1  /i t  (closed  loop  system) 

Created  by  OCMFRKQl  and  OC.MFRFQJ  (data  files): 

yp  —  human  operator  ye  =  controlled  system 

upyc  =  loop  transfer  function  p/»i  —  remnant 


3.2  Using  Program  CC 

Program  CC  is  a  command  driven  computer  aided  -  control-system-design 
package.  It  is  written  in  compiled  BASIC  and  operates  under  the  DOS 
operating  system  on  the  IBM- PC  and  compatible  personal  computers.  Pro¬ 
gram  CC  works  with  linear  systems,  either  analog  or  digital,  which  are 
modeled  with  either  transfer  functions  or  state  space  equations.  A  large 
number  of  classical,  sampled  data,  and  state  space  algorithms  are  imple¬ 
mented,  including  an  extensive  amounts  of  interactive  graphics.  The  user 
environment  in  Program  CC  is  robust,  friendly,  and  very  powerful. 

Program  CC  uses  macron  and  user -defined  commands  to  tailor  its  use 
to  a  particular  problem.  The  SETIP  command,  as  explained  here,  is  used 
to  make  sure  that  the  program  has  access  to  these  features. 

Version  4  of  Program  CC  must  be  used  for  the  optimal  control  model 
of  the  human  operator.  An  introduction  to  using  Program  C'C,  Version  4  is 
contained  in  Appendix  A.  This  introduction  concentrates  on  the  parts  of  the 
program  used  for  the  optimal  control  model.  A  more  extensive  introduction 
is  contained  in  the  Tutorial  Manual  [3] 

3.2.1  Setup 

Obtain  a  copy  of  Program  CC.  Version  4  and  then  load  it  onto  the  hard 
disk  of  an  IBM  PC  or  compatible  personal  computer.  The  program  con¬ 
sists  of  about  40  separate  modules  with  the  names  CC.EXE.  CCT.EXE.  and 
so  on.  Place  all  of  these  modules  into  a  subdirectory  named  CC.  Create  a 
subdirectory  named  CC’/DATA  to  store  the  data,  and  a  subdirectory  named 
CC/OC.M  to  store  the  OCM  macros.  Copy  the  OCM  macros  listed  in  Ap¬ 
pendix  1!  into  the  CC/OC.M  directory.  The  OCM  macros  work  only  with 
Program  CC.  Version  4.  and  will  not  run  properly  with  Version  3. 

Program  CC  has  the  ability  to  locate  data,  macros,  and  modules  on 
different  drives.  Because  of  the  large  number  of  modules  (i.e.  overlays)  used 
by  Program  CC  the  efficiency  of  operation  is  improved  by  placing  some  ol 
these  in  a  RAM  disk.  The  following  modules  which  are  used  by  the  OCM 
are  listed  in  order  of  priority: 
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CC 

root  program 

CC1 

macro  processor  and  utilities 

CC25 

state  space  equations 

CC31 

eigenvalue  calculations 

CC'34 

Lyapunov  solver 

CC2 

transfer  function  display- 

CC3 

polynomial  root  finder 

CC4 

transfer  function  equations 

All  of  the  above  modules  will  fit  into  a  RAM  disk  of  640  Kbytes.  It  is 
also  possible  to  put  the  data  used  for  the  OC’M  into  the  same  RAM  disk, 
in  which  case  it  should  bo  expanded  by  at  least  128  Kbvtes.  For  safety,  in 
case  of  a  power  failure  or  a  system  crash,  it  is  better  to  place  the  data  in  a 
subdirectory  on  the  hard  disk.  Data  stored  on  the  hard  disk  can  is  not  lost 
when  Program  CC  is  stopped. 

Use  the  CC/SETUP  command  to  establish  the  pointers  to  the  above 
items  will  be  demonstrated  after  a  few  more  preliminary  comments.  How 
to  setup  a  RAM  disk  depends  on  your  particular  hardware  configuration 
and  version  of  DOS.  Commands  similar  to  the  following  should  be  placed 
in  your  CONFIG.SYS  file  (the  files  and  buffer  commands  are  not  related  to 
the  RAM  disk  but  are  nevertheless  a  good  idea): 

device=c:\dos\edisk.sys  610  128  512/e 

files =20 

buffers=32 

3.2.2  Executing  Macros 

A  macro  is  an  indirect  file  containing  Program  CC  commands.  The  com¬ 
mands  arc  listed  in  a  macro  exactly  as  they  would  be  if  entered  from  the 
keyboard.  Macros  can  be  nested  and  parameters  can  be  inserted.  Macros 
are  executed  by  including  the  (*  symbol  before  t he  macro  name,  for  example: 

CC>  QQCMLQR,  P40,  1,  le-3,  5 

The  macros  used  for  t he  human  optimal  control  model  are  all  written  so 
that  they  echo  the  parameters  and  then  pause  for  a  response.  Press  the  FI 
function  key  to  abort,  and  any  other  key  to  continue.  The  Fl  function  key 
can  be  pressed  at  any  time  to  stop  execution  (though  it  is  not  recognized 
during  an  overlay  swap  and  sometimes  several  key  strokes  are  required). 
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3.2.3  U ser  Deft ned-Coin mauds 


A  user -dcjintd --command  is  a  BASI<  program  which  is  connected 

to  Program  CC  during  oxeru* ion .  On-  Mich  user  defined-command  is  used 
here,  OCM.EXE.  which  has  2  mam  functions: 

1.  Controls  the  I.QR  iterations  by  automatically  changing  the 
input  weight  g. 

2.  Controls  the  KBF  iterations  by  automatically  updating  the 
observation  and  motor  noise  intensities.  Also  computes  the 
erfc  function  by  means  of  a  piecewise  linear  approximation, 
and  prints  a  summary  at  the  end  of  each  KBF  iteration. 

3.2.4  Example 

Start  the  program  from  the  DOS  level  by  typing  CC.  Execute  the  SETUP 
command  as  shown  below  to  establish  pointers  to  the  data,  macros,  and 
modules  used  by  the  OC'.M: 

CC>setup  SETUP  command 

2  establish  pointers  to  data,  macros,  and  pro¬ 

gram  modules 

2  use  automatic  option 

\cc\data  data  location 

\cc\ocra  macro  location 

\cc\data  more  data  (not  used  by  the  OCM) 

d:  name  of  RAM  disk 

4  establish  a  link  to  the  user  defined  command 

1  add  a  command 

ocm  name  of  command 

ocm  name  of  BASIC  program 

n  no  more  commands 

4  return  to  main  menu 

5  change  the  SSSETUP  file 
CC>EXIT  leave  the  program 

Make  a  back  up  copy  named  SETUP  of  the  file  just  created.  After 
completing  the  initialization  start  Program  C'C  with  a  DOS  batch  file  such 
as  the  following,  which  copies  the  SETUP  file  and  the  first  three  of  the 
recommended  modules  into  the  D:  RAM  disk: 
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COPY  SETUP  $$ SETUP 
COPY  CC .EXE  D: 

COPY  CC1.EXE  D: 

CUPY  CC25.EXE  D: 

CC 

3.3  The  Controlled  System 

Ihe  controlled  system,  denoted  as  ( Yc)ati •  includes  the  controlled  element 
Yc(s),  the  noise  dynamics  Vu,(s),  and  the  interconnection  structure.  The 
inputs  of  {Yc)aii  are  the  operator  input  u  and  t lie  external  noise  source  w. 
The  outputs  are  the  error  e  and  the  error  rate  c.  The  structure  depends  on 
where  the  noise  is  injected,  with  the  three  possible  choices  (for  a  single-input 
system)  being  at  the  output  of  Fc(s),  the  input,  or  somewhere  between.  Use 
the  macros  OCMYC1.  OC'M^  C'2,  and  OCMYC'3  to  create  a  state  space 
quadruple  for  {Yc)au  starting  from  transfer  function  descriptions  of  Yc(s) 
and  Yw(s).  The  block  diagrams  in  Figure  4  explain  the  three  cases. 

3.4  The  Optimal  Solution 

After  creating  the  model  of  the  controlled  system  the  I.QR  and  KBF  optimal 
control  problems  must  be  solved.  Fach  of  these  is  an  iterative  process.  The 
0(  MALL  macro  defines  the  following  default  parameter  values  and  then 
controls  both  iterations: 

q  =  1 
r,v  =  .1 

Tthrrsh  —  .001 

Vu-  =  1 

T  =  .2 

Py  i  —  Py 2  *  Pua  =  —  2 0(1  B 
7 1  =  72  =  0 
/=  1 

(l  Bt htrr  sh  —  -1 

1  he  OCMALL  macro  works  by  calling  several  lower  level  macros,  two  of 
which  call  the  OCM  user  defined  command  (lTDC): 
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OCMLQR 

OCMSETUP 

OC’MKBF 

OCMPILOT 


solves  the  F.QR  problem,  calls  OC'M 
enters  parameters  used  by  OCMKBF 
solves  the  KHF  problem,  call  OCM 
computes  Vp( s ) 


The  default  parameters  listed  above  can  be  changed  by  adjusting  the  inputs 
to  the  lower  level  macros.  If  these  parameters  are  changed  then  change  the 
name  of  OCMALL  to  something  like  OCMALL1. 

The  OCMLQR  macro  changes  the  input  weight  g  and  repeats  the  LQR 
problem  until  the  achieved  r /v  is  within  Tthrtsh  of  the  desired  value.  The 
following  ad-hoc  method  is  used:  ( 1 )  start  with  an  initial  guess.  (2)  increase 
or  decrease  by  a  factor  of  10  until  the  desired  r.y  is  surrounded,  (3)  con¬ 
verge  using  a  binary  search.  After  convergence  you  arc  presented  with  the 
following  choices: 


1-stop 

2=1  more  iteration 
3=change  threshold 

Use  the  first  option  to  stop  the  LQR  iteration  and  continue  to  the  KBF 
problem.  Use  the  second  option  to  continue  the  LQR  iteration  one  more 
step,  and  use  the  third  option  to  set  a  more  precise  threshold.  The  LQR 
iterations  typically  go  very  fast. 

The  OCMSETUP  macro  sets  up  the  KBF  iterations  by  storing  the  re¬ 
quired  parameters  in  set  locations.  The  OCMKBF  macro  controls  the  KBF 
iterations,  which  is  a  longer  and  more  difficult  problem  then  the  LQR  itera¬ 
tions.  The  KBF  iterations  work  bv  changing  V'vi .  \'y2.  and  VUa,  as  described 
in  the  technical  write-up.  until  the  desired  noise  ratios  are  achieved. 

At  the  end  of  each  KBF  iteration  a  summary  is  presented  of  the  time  do¬ 
main  performance  measures:  optimal  cost,  the  mean  square  errors,  the  noise 
intensities,  and  the  dB  noise  ratios  (see  the  examples).  Any  of  these  can  in 
theory  be  used  as  a  stopping  criteria,  however  the  one  which  is  implemented 
is  for  all  of  the  noise  ratios  to  be  within  dBthrtsh  of  their  desired  values.  You 
can  judge  the  progress  of  the  iteration  from  the  last  line  of  the  summary, 
which  give  the  current  value  of  the  maximum  noise  ratio  difference. 

A  different  threshold  may  be  used  depending  on  the  computation  time, 
the  rate  of  convergence,  and  your  patience:  but  more  than  1  dB  is  not  rec¬ 
ommended.  Make  the  change  by  adjusting  the  c IBtkrtsh  parameter  in  the 
OCMSETUP  macro.  The  number  of  required  iterations  depends  on  the  ini¬ 
tial  noise  intensity  estimates,  the  dynamics  of  the  system  (hard  to  control 
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systems  and/or  unstable  systems  tend  to  require  more  iterations),  tlio  atten- 
tional  fraction  and  indifference  threshold  levels,  and  the  stopping  criteria. 
Anywhere  from  3  to  more  than  50  iterations  may  be  needed.  Typically  less 
than  5  iterations  are  required  to  reach  a  noise  ratio  window  of  1  dB,  but 
convergence  from  that  point  on  can  be  slow.  The  best  way  to  significantly 
reduce  the  number  of  iterations  is  to  give  good  starting  values  for  Vy\ ,  Vy2, 
and  V'Uo.  As  a  practical  matter,  this  can  only  possible  if  the  same  problem  or 
a  similar  problem  has  been  previously  solved.  If  convergence  is  not  achieved, 
which  sometimes  happens,  then  stop  the  process  using  the  Fl  function  key. 
Even  after  stopping  it  is  still  possible  to  continue  with  further  analysis. 

When  the  KBF  iteration  is  completed  you  are  presented  with  the  follow¬ 
ing  choices: 

l=stop 

2=1  more  iteration 
3=change  threshold 
4=line  printer  listing 

The  last  option  is  strongly  recommended,  which  creates  a  line  printer  list¬ 
ing  of  the  last  iteration  summary.  Permanent  records  are  always  desirable, 
though  in  the  heat  of  the  moment  sometimes  forgotten.  Use  the  first  op¬ 
tion  to  st op  the  KBF  iterations  and  continue  with  the  analysis.  The  second 
option  continues  with  one  more  step,  and  the  third  option  prompts  for  a 
smaller  t hreshold. 

The  OCMPILOT  macro  creates  a  state  space  approximation  of  Up(s). 
The  information  needed  to  compute  Yp(s)  is  ( 1 )  the  controlled  system  (Yc)alh 
(2)  the  LQR  solution,  decomposed  into  the  Tv  and  Lm  matrices.  (3)  the  KBF 
gains  H],  and  (4)  the  visual  delay  r.  The  state  space  model  created  by  the 
OCMPILOT  macro  is  an  approximation  because  a  Pade  approximation  is 
used  in  place  of  an  exact  delay.  A  second  order  approximation  is  recom¬ 
mended.  The  parameters  for  the  OCMPILOT  macro  are  the  delay  and  the 
order  of  the  approximation.  The  remainder  of  the  information  is  assumed 
to  be  in  the  locations  used  by  the  previous  macros.  The  default  name  for 
the  state  space  model  of  Yp(s)  is  P25. 

3.5  Analyzing  the  Results 

After  executing  the  OCMALL  macro,  or  its  constituent  parts,  there  are  sev¬ 
eral  possible  directions  to  proceed  with  the  analysis.  The  suggestions  here 
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are  not  meant  to  be  exhaustive.  Think  of  the  macros  mentioned  here:  0C- 
MALL,  OCMG,  OCMFREQl.  and  0CMFREQ2,  as  objects  to  be  changed 
according  to  the  specific  task  at  hand.  It  is  of  course  best  to  change  the 
names  if  modifications  are  made. 

The  time  domain  performance  measures  are  listed  at  the  end  of  each 
IvBF  iteration.  The  parameters  are  stored  in  matrices  as  listed  in  Tables  3. 
These  parameters  may  be  the  end  point  of  your  analysis,  in  which  case  no 
further  work  is  needed.  For  example: 

1.  The  optimal  cost  J  can  be  used  to  obtain  a  pilot  opinion  rating.  The 
value  of  J  is  listed  in  the  OCM  summary.  By  referring  to  Table  3  it 
is  seen  that  J  is  also  stored  in  the  P23  matrix.  The  definition  of  J  is 
listed  in  Table  2;  if  additional  normalization  is  needed  multiply  J  by 
the  appropriate  scale  factor. 

2.  The  optimal  mean  square  tracking  error  crj^  can  be  compared  against 
experimental  results.  The  value  is  listed  in  the  OCM  summary,  along 
with  other  mean  square  results,  and  is  also  stored  in  the  P20  matrix. 

The  OCMFREQl  macros  use  the  state  space  approximations  for  Vp(s) 
and  y'c,(s)  to  compute  frequency  data  files  for: 


Fp(-M 


<l>nn,(s)l/2 

The  macro  ends  with  a  Bode  plot  created  with  the  BLOT  command.  The 
square  root  of  the  remnant  is  computed,  as  is  customary.  (The  dB  scale 
used  in  Program  CC  is  201og10.  Plotting  the  square  root  of  the  remnant 
is  the  same  as  using  a  power  dB  scale  of  101og]0.)  The  cursor  can  be  used 
to  determine  frequency  domain  parameters  such  as  bandwidth  and  phase 
margin. 

The  0CMFREQ2  macros  differ  in  that  they  use  an  exact  calculation  of 
e_ST  as  an  intermediate  step  for  the  same  frequency  data  files  as  above.  I 
have  not  yet  come  across  an  example  where  there  is  any  significant  difference 
in  the  crossover  region  between  the  two  methods. 

Another  type  of  analysis  is  to  obtain  a  transfer  function  approximation 
of  Tp(  s),  which  can  then  be  used  to  obtain  low  order  approximations.  Clas¬ 
sical  (i.e.  structural)  pilot  model  parameters  such  as  pilot  lead  and  effective 
delay  are  easiest  to  obtain  using  low  order  models  of  Vp(s).  After  using  the 
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OCMPILOT  macro  to  obtain  a  state  space  model  for  Yp(s),  use  the  OCMG 
macro  to  convert  to  a  transfer  function  model.  The  optimal  control  model  of 
the  human  operator  includes  the  negative  feedback  sign  with  Y'p(s),  which  is 
the  opposite  convention  used  by  the  classical  human  operator  model.  Mul¬ 
tiply  the  >'p(5)  computed  by  the  OCMG  by  -1  to  switch  to  the  convention 
used  by  the  classical  model.  The  OCMG  macro  also  computes  I’r(.s),  in  case 
it  is  not  already  available.  Multiply  Yp(s)  and  Vr(s)  together  to  obtain  the 
loop  transfer  function,  which  can  then  be  checked  against  the  K /s  crossover 
rule. 

The  NEAR  command  is  used  in  the  OCMG  macro  to  cancel  poles  and 
zeros  of  Yp($)  and  lrc(s).  There  will  be  considerable  cancellation,  and  more 
may  be  desired.  For  further  cancellation  use  the  NEAR  command  with  larger 
tolerances.  The  LFAPPROX  (low  frequency  approximation)  command  can 
be  used  to  truncate  poles  and  zeros  larger  than  a  threshold  frequency,  with 
the  option  of  replacing  the  effective  delay  of  the  truncated  modes  with  a 
Pade  approximation  of  the  delay. 

Converting  to  transfer  functions  is  probably  the  best  way  to  quickly  ana¬ 
lyze  the  results  of  the  OCM,  because  Program  CC’s  user  interaction  is  faster 
and  more  convenient  in  the  transfer  function  domain.  In  particular,  the  plot 
options  can  be  used  for  precise  bandwidth  and  phase  margin  calculations. 

3.6  Computation  Time 

Three  factors  are  important: 

Number  of  states  in  Yc 
Number  of  KBK  iterations 
Micro- processor  clock  speed 

The  following  operations  in  each  KB1  iteration  take  the  most  time,  each  of 
which  is  a  n3  operation,  where  »  is  the  number  of  states  in  Yc: 

Riccati  equation  solution 
mat  rix  exponential 
Lyapunov  equation  solution 

It  is  best  to  use  low  order  closed  loop  effective  models  for  the  controlled 
system  Tc,  and  to  use  only  1st  and  2nd  order  filters  for  the  driving  noise.  If 
Y  has  5  or  6  states,  then  on  an  AT  with  a  8  MHz  clock  it  will  take  about 
50  seconds  for  each  KBF  iteration.  Due  to  the  n  *  dependence,  doubling  the 
number  of  states  requires  8  times  as  many  computations. 
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The  number  of  inputs  is  not  the  dominant  factor  in  computation  time 
(probably  a  m2  dependence,  where  m  is  the  number  of  inputs).  Systems 
with  two  or  three  inputs,  however,  typically  have  two  or  three  times  as 
many  states,  especially  if  the  systems  are  independent  in  each  axis. 

The  computation  time  is  proportional  to  the  number  of  KBF  iterations, 
obviously.  How  many  iterations  are  required  is  discussed  in  Section  3.2.4 
Computation  time  is  just  as  obviously  proportional  to  the  micro-processor 
clock  speed.  Use  your  fastest  computer. 

Overlay  and  data  read  times  can  be  very  significant  if  a  HAM  disk  is  not 
used.  It  is  definitely  worth  the  trouble  learning  how  to  use  a  RAM  disk. 


4  Examples 

Duplicate  the  following  (X  'M  examples  to  gain  experience. 

4.1  Integrator 

The  controlled  system  is  a  simple  integrator  in  this  example  taken  from  [1]. 
The  input  of  the  system  is  driven  bv  colored  noise,  the  human  operator  visu¬ 
ally  determines  the  error  and  error  rate,  and  manually  controls  the  input  so 
as  to  minimize  the  mean  square  error.  The  example  has  following  controlled 
system  and  parameters: 


Driving  noise:  V'u,  =  8.8.  resulting  in  E{rf}  =  2.2 
LQIt :  <7  =  1.  g  =  .00017.  resulting  in  t\  =  .OS. 

Visual  delay:  r  =  .15 

Noise  ratios:  pVl  =  .01,  =  .01.  =  .00516  (-20.-20,-25  d B ) 

Indifference  thresholds:  t\  =  0.  tj  =  0 
Fractional  attention:  /  =  1 
Order  of  \ '  :  n  =  2 

I  he  parameters  r.v.  r.  and  pUa  were  adjusted  by  the  authors  [1]  in  or¬ 
der  to  obtain  a  close  match  with  experimental  results.  The  Program  f’C 
implementation  agrees  with  all  of  the  published  results  except  for  the  rem¬ 
nant  frequency  response  (its  not  clear  whether  or  not  the  remnant  formulas 
match ). 

The  state  space  description  of  (>’  }.,//  is  provided  in  the  problem  state 
ment.  and  therefore  can  be  directly  entered: 

STATE>p40=( -2,0,0, 1 ;  1.0, 1,0;  0,1, 0,0;  1,0, 1,0) 
STATE>p40=CHST (p40 ,2) 

Ihe  CHS  I  function  is  used  to  convert  a  real  matrix  into  a  state  space 
quadruple,  in  this  case  with  2  states.  It  is  helpful  to  compute  the  same 
result  starting  from  transfer  functions: 

CC>yc=l/s 
CC>yw=l/ (s+2) 


COCOCMYC2,yc  ,yw  ,p40 
CC>p40*p40(s ,(2,1)) 

In  this  example  the  noise  in  injected  in  the  input  of  the  plant,  hence  the 
OCMYC2  macro  is  used  to  create  P40.  The  odd  looking  last  command 
switches  the  order  of  the  2  states,  an  optional  step  which  is  used  to  obtain 
the  same  state  space  realization  as  given  in  the  problem  statement. 

Now  compute  the  LQR  problem.  The  value  of  g  is  already  provided,  so 
that  no  iteration  is  required: 

CC>CLQR,p40,l, .00017, .08,. 001 

The  rv  actually  achieved  is  .08074,  which  falls  within  the  threshold  of  .001. 

There  are  many  parameters  required  for  the  KBF  problem.  List  the 
parameters  by  running  the  OC’MSETl'P  macro  followed  by  the  FI  function 
key,  and  then  include  the  numbers  on  the  command  line.  Carefully  check 
before  proceeding  to  the  KBF  iterations: 

CC>eOCMSETUP 

FI 

CO0OCMSETUP,  8.8,  0.15,  -20,  -20,  -25,  1,1, 1,0, 0,1 
CO0OCMKBF,  .1 

Six  iterations  are  required  to  achieve  the  desired  noise  ratios.  One  more 
iteration  is  computed  for  good  measure: 

Table  1:  KBF  Iterations 


m 

1 

-16.23 

2 

13.47 

-15.33 

3 

17.23 

-18.72 

-24.32 

4 

-19.23 

-19.68 

-24.84 

5 

-19.81 

-19.92 

-24.96 

6 

-19.96 

-19.98 

-24.99 

7 

-19.99 

-20.00 

-25.00 

The  iteration  summary  after  the  7th  iteration  is  shown  in  Figure  5.  The 
following  noise  intensities  which  were  converged  to: 

\V2  =  .09087  rUo  =  .  01815 


Vyl  =  .00371 


Use t  his  hard  achieved  informal  ion  to  speed  the  solution  of  t  ho  problem  the 
next  time  around: 


COQOCMSETUP,  8.8,  0.15,  -20,  -20,  -25,  .00371, 

.09687,  .04815  ,0,0,1 
COQOCMKBF,  .1 

Anytime  you  want  to  bring  the  iteration  summary  back  onto  the  screen 
simply  call  the  OCMKBF  macro  once  more.  Begin  the  analysis  by  noting 
the  time  domain  performance  parameters  listed  in  the  summary  of  the  last 
iteration: 


UV2  “ 

<  =  3.86 
J  =  .159 

Run  the  OCMI’ILOT  macro  to  obtain  a  state  space  model  for  Yp(s). 
Follow  this  with  the  GEP  command  (GEP  stands  for  the  Generalized  Eigen¬ 
value  Problem),  which  creates  a  9th  order  transfer  function  model  for  }p(s): 

COQOCMPILOT ,  0.15,  2 

COSTATE 

STATE>GEP ,p24 ,yp 

Change  the  sign  of  >p(.s)  to  conform  with  the  classical  convention.  Use 
the  NEAR  command  to  cancel  poles  and  zeros  with  an  absolute  difference 
less  than  10-4,  and  then  use  the  LFAPPROX  command  to  replace  all  of  the 
dynamics  greater  than  5  rad/sec  with  a  1st  order  delay  approximation: 

STATE>cc 
COyp— yp 

CC>NEAR,yp,ypl , 1 , le-4 
CC> LFAPPROX ,ypl ,yp2 ,5,1 

Figure  5  shows  the  V'.(s)  and  )  , (.(.s)  transfer  functions,  as  well  as  the  transfer 
functions  for  V'p(.s)  before  and  after  cancellation.  The  final  low  frequency 
approximation  for  V'p(.s)  is: 


Iteration  #  7 

Optimal  cost:  q,  g,  output,  input  rate,  total 

1 . 0000E+00  1 . 7000E-04  1.1803E-01  4.1207E-02  1.5923E-01 

Performance  :  E{y_lA2),  E(y_2'2),  E{u_a'2),  E(uA2),  E((du/dt)A2) 

1 . 1803E-01  3 . 08  34E+00  4.8469E+00  3.8633E+00  2.4240E+02 

Old  noise  intensities  (V_yl,  V_y2,  V_ua) :  3.7167E-03  9.6967E-02  4.8176E-02 

New  noise  intensities  (V_yl,  V_y2,  V_ua) :  3.7079E-03  9.6869E-02  4.8152E-C2 

Noise  ratios  dB  (rho_yl,  rho_y2,  rho_ua) :  -19.9897  -19.9956  -24.9978 

Max  noise  ratio  difference  =  1.028061E-02  dB,  Threshold  =  .1  dB 


Controlled  element  and  noise  models 
1 

YC(s)  =  — 


1 

YW ( s )  =  - 

s+2 


Pilot  model  and  low  order  approximations 

179 . 1459 (  0)  (  2)  (  3.252133)  (  6.386982)  (  12.38519)  (  12.56791) 
[-.8660254,  23.09401] 

YP(s)  - - - - - - 

(  0)  (  1.990881)  (  2)  (  6.460116)  (  12.37807)  (  12.38519) 

[  .3671649,  23 . 26644 ] (  42.51991) 


17  9 . 14  59  (  3.252133)  (  6.386982)(  1 2 . 5 679 1 )  [ - . 86602 54 ,  23.09401] 

YPl(s)  =  - 

(  1.990881)  (  6.460116)  (  12.37807)[  .  3671649,  23.26644  ] 

(  42.51991) 


-4 . 16696  (  3.252133)  (-15.44067) 

YP2 (s)  =  - 

(  1.990881) (  15.44067) 


Figure  5:  Summary  of  (he  IntegraJ  OCM  Example 


This  is  similar  to  what  a  classical  (i.e.  structural)  pilot  model  would  predict, 
except  that  the  optimal  V'p(s)  contains  the  extra  lag  term  (s  +  3.3)/(s  +  2.0). 
TI,p  effect  of  this  tprm  is  to  increase  the  low  frequency  gain,  which  improves 
the  tracking  response  to  low  frequency  inputs. 

A  frequency  plot  of  the  important  transfer  functions  is  obtained  by  the 
following  macro  call: 

COSOCMFREQ1,  .01,100,100 

The  resulting  Bode  plot  shown  in  Figure  6.  The  parameters  specify  100 
points  from  .01  to  100  rad/sec.  The  phase  of  Yp(s)  and  }pVc(.s)  is  180°  away 
from  that  expected  from  the  classical  pilot  model  convention.  Several  of  the 
standard  frequency  domain  performance  parameters  are  listed  below.  (Use 
the  MARGIN  and  POINT  commands  to  help  with  some  of  the  following 
numbers). 


-A- 

OPM 

Te 


4.88  rad /sec 
41.9° 

it  7 r 

. 1 7  sec 


The  commands  used  to  duplicate  this  example  are  stored  in  the 
EX  1. MAC  listed  in  Appendix  B. 


40 


V_p,  V_c .  V_p»il_c ,  and  Pht 


10  2  10  1  I00  10 1  10Z 


Frequency  (rad/sec) 


Figure  6:  Frequency  Responses  for  Integral  OCM  Example 
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4.2  Double  Integrator 

Another  standard  human  operator  tracking  problem  uses  a  double  integra¬ 
tor.  In  this  example  inject  the  driving  noise  at  the  output  of  the  controlled 
element,  and  model  the  noise  as  unit  intensity  white  noise  passed  through  a 
2nd  order  Butterworth  filter  with  a  break  at  .5  rad/sec.  Use  defaults  for  all 
of  the  remaining  OCM  parameters: 

Controlled  element:  Vc(s)  =  1/s2 
Driving  noise:  Yw(s)  =  l/(4.s2  +  2.83  *  s  +  1) 

Driving  noise:  Vu.  =  1 
LQR:  q  =  1,  t\  =  .1. 

Visual  delay:  r  =  .2 

Noise  ratios:  =  .01.  pyi  -  .01./>„o  =  .01  (-20.-20,-20  dB) 

Indifference  thresholds:  t\  -  0.  0  =  0 
Fractional  attention:  /  =  1 


Solve  the  problem  with  the  following  commands: 

CC> yc=l/s'2 
COBUTTER, yw,  .5,2 
COQOCMYCl  ,yc  ,yw,p40 
CC>QOCMALL,p40 

Obtain  transfer  function  approximations  of  Vp(.s)  with  the  following  ad¬ 
ditional  commands: 

STATE>GEP ,p24 ,yp 
STATE>CC 

CC>NEAR,yp,ypl , 1 , le-4 
COLFAPPROX ,ypl  ,yp2,5, 1 
CONEAR, yp2, yp3, 1 ,  .2 

The  final  result  is: 

YPU »)  =  2.1(s+  1 .3  )o~iia 

The  iteration  summary  after  16  iterations  is  shown  in  Figure  7,  together 
with  the  transfer  function  results.  The  frequency  responses  are  shown  in 
Figure  8.  The  transfer  functions  were  used  to  obtain  this  plot,  and  hence  the 
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straightline  asymptotes  can  he  included.  The  frequency  domain  performance 
parameters  are  listed  below: 

ojc  =  3.16rad/sec 
0PM  =  26.3° 
r-  = 

=  .35  sec 

The  commands  used  to  duplicate  this  example  are  stored  in  EX2.MAC. 
The  LQR  and  KBF  problems  are  started  at  the  parameters  which  have 
already  been  converged  to. 


Iteration  #  16 

Optimal  cost:  g,  g,  output,  input  rate,  total 

1 . OOOOE+OO  6 . 4053E-05  1.6305E-02  8.0904E-03  2.4395E-02 

Performance  :  E{y_lA2),  E(y_2A2),  E(u_a*2),  E(u*‘2),  E(  (du/dt)  "2 ) 

1.6305E-02  1 . 1800E-01  2.6566E+00  2.2338E+00  1.2631E+02 

Old  noise  intensities  (V_yl,  V_y2 ,  V_ua) :  5.1638E-04  3.7338E-03  8.4O74E-02 

New  noise  intensities  (V_yl,  V_y2,  V_ua) :  5.1223E-04  3.7070E-03  8.3460E-02 

Noise  ratios  dB  (rho_yl,  rho_y2,  rho_ua):  -19.9650  -19.9688  -19.9682 

Max  noise  ratio  difference  =  3.499794E-02  dB,  Threshold  *  .1  dB 


Controlled  element  and  noise  models 
1 

YC(s)  =  - 

s*2 


1 

YW(s)  =  - 

4sA  2  +2.83S+1 


Pilot  models  and  low  order  approximations 

500 (  0) *2  (  .714,  . 438 ] [  .707,  .5](  1.33) (  3.11)(  10) (  10.1) 
[-.866,  17.3] 

YP(S)  . . . - . — . - . 

(  0) "2  [  .707,  . 5 ] * 2  {  3 . 1)  (  9.99)(  10) [  .275,  10.6)[  .865,  22.4] 


500 [  .714,  . 438 ]  (  1.33)(  3.11)  (  10.1)[-.866,  17.3] 

YPl(s)  =  - 

[  .707,  . 5 ] {  3 . 1 ) (  9.99) [  .275,  10.6][  .865,  22.4] 


-2 . 68  [  .714,  .438]  (  1.33)  (  3.11)  (-8.69) 

YP2 ( s )  =  - 

[  .707,  .5] (  3.1) (  8.69) 


-2 . 05 (  1.33) (-8.69) 

YP3 ( s )  =  - 

(  8.69) 


Figure  7:  Summary  of  Double  Integrator  OCM  Example 


10  1  10a  10 1  10Z 


Frequency  (rad/sec) 


Figure  S:  Frequency  Responses  for  Double  Integrator  OCM  Example 
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4.3  Single  Axis  Dander  Example 

Dander  [15]  has  reported  experimental  findings  from  single  and  multi-axis 
tracking  tasks.  Dynamically  independent  single,  two-,  and  three-axis  track¬ 
ing  experiments  were  conducted,  and  then  subjective  pilot  opinion  ratings 
(PORs)  were  given  for  each  task  using  the  Cooper-Harper  rating  procedure. 
(This  procedure  yields  a  POR  ranging  from  1  (best)  to  9  or  10  (worst), 
which  is  dominated  by  the  pilot's  mental  workload  required  to  achieve  the 
performance  implied  by  a  given  mission  phase.)  One  of  the  objectives  of 
the  original  experiment  was  to  predict  multi-axis  PORs  based  on  single  axis 
results.  The  best  way  to  make  this  extension  is  still  a  open  question,  and 
the  Dander  data  is  still  the  best  data  base  upon  which  to  test  out  theories. 

The  first  part  of  this  example  solves  the  OCM  for  a  single  axis  task,  and 
the  second  part  uses  the  OCM  to  predict  single  and  multi-axis  PORs.  The 
method  used  to  determine  the  optimal  multi-axis  cost  from  the  single  axis 
optimal  costs  is  described  in  Volume  2  of  this  report.  More  explanation  of 
the  Dander  data  and  POR  predictions  is  provided  by  McRuer  and  Schmidt 
[15]  and  Volume  2  of  this  report.  Further  background  on  the  use  of  optimal 
cost  for  PORs  in  contained  in  [13]. 

The  single  axis  controlled  system  and  driving  noise  filter  are  listed  below: 

-l(-(>  1)(  -9) 

(0)[.7..25](5) 

.2219 
[.7.  .5 

The  driving  noise  is  added  to  the  output  of  Tc(.<).  The  mean  square  output 
error  due  to  the  driving  noise  is  =  .14.  The  experimental  parameters  are 
listed  below,  note  that  nonzero  indifference  thresholds  are  used: 

f  XL  7~V  {)  y  |  f)y  2  Pua _ (j _ 0 _ f 

1  .2  .1  .01  .01  .01  .015  .025  1 

Enter  the  transfer  functions,  combine  them  to  form  the  controlled  system, 
and  then  solve  the  OCM  problem: 

COG0CMYC1 ,  yc  ,yw  ,p40 
CC><JOCMALL,p40 

I  he  following  time  and  frequency  domain  performance  parameters  have  been 
computed  using  the  OCM  solution: 


T,(.s)  = 

ru.(s)  = 


40 


.00018  .005  .049  .16  .20 


si _ ■!_ 

9.8  .0067 


Of'M  Tt 

3.2  r/s  37°  .29  sec 

See  the  EX3  macro  in  Appendix  B  for  more  of  the  commands  which 
duplicate  this  example.  The  converged  to  values  of  g ,  Vyl,  Vy2,  and  VUa  are 
used  in  this  macro,  so  that  only  1  iteration  is  needed  for  each  of  the  LQR 
and  KBF  problems.  Figure  9  shows  the  iteration  summary  and  several 
difference  approximations  of  Vp(s).  Figure  10  shows  a  Bode  plot  of  the 
important  transfer  functions,  complete  with  straightline  asymptotes.  The 
low  order  approximation  for  the  pilot  model  is: 

=  2.6(5+  1.4)(5  +  5.2)  21a 

p  (s  +  .9) 

The  classical  pilot  model  is  Yp  =  1  .l(.«+5)e_-25j.  The  lead  is  used  to 
maintain  a  K/s  like  crossover,  and  the  delay  is  determined  as  a  function 
of  the  3  r/s  crossover.  The  OCM  derived  Y'p,  despite  its  high  order  (n+5, 
where  n  is  the  order  of  Y'c  plus  the  order  Yw),  is  basically  the  same  around 
crossover,  but  differs  at  lower  frequencies  by  including  trim  terms,  and  differs 
at  higher  frequencies  by  including  an  ill  understood  collection  of  terms  which 
can  effectively  be  grouped  into  a  delay.  The  much  vaunted  nemo  muscular 
mode  is  cancelled  by  the  OCM  solution  and  does  not  appear  in  V'p.  It  is 
nevertheless  important  to  set  the  neuro-muscular  mode  as  part  of  the  LQR 
solution,  because  this  method  of  selecting  the  quadratic  weight  determines 
the  closed  loop  bandwidth. 


4.4  Multi-Axis  Dander  Example 


The  Dander  experiment  varied  the  dynamics  in  each  of  3  axes  to  create  a 
large  number  of  combinations.  Here  only  a  single  set  is  used:  6h<Pl3h  in 
Dander’s  notation.  The  OCM  is  solved  for  each  axis,  and  then  the  optimal 
costs  for  each  axis  are  combined  as  demonstrated  below  to  obtain  an  optimal 
cost  for  the  multi-axis  case.  After  all  of  these  optimal  costs  are  determined 
then  predictions  of  PORs  are  made. 

The  0^  axis  (axis  #1)  was  considered  in  the  previous  subsection.  The 
full  attention  case  for  the  <i>i  axis  (axis  #2)  is  summarized: 


Ye(s) 


0-5(0.1) 

( 1.5 )[— .84,0.5] 


13.3 

[,7,.5] 
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Iteration  #  1 

Optimal  cost:  q,  g,  output,  input  rate,  total 

1 . 0000E+00  1 . 8  000E-04  4.9850E-03  1.7647E-03  6.7497E-03 

Perrormance  :  E{y_l'2),  E{y_2'2),  E{u_a'2),  E(u*2),  E({du/dt)'2) 

4 . 9850E-03  4.8742E-02  1.9835E-01  1.6504E-01  9.8040E+00 

Old  noise  intensities  (V_yl,  V_y2,  V_ua) :  2.7070E-04  2.0260E-03  6.2760E-03 

New  noise  intensities  (V_yl,  V_y2 ,  V_ua) :  2.6840E-04  2.0101E-03  6.2313E-03 

Noise  ratios  dB  (rho_yl,  rho_y2 ,  rho_ua) :  -19.9630  -19.9657  -19.9689 

Max  noise  ratio  difference  =  3.701782E-02  dB,  Threshold  =  . 1  dB 


Controlled  element  and  noise  models 
4<  -04) (  .9) 

YC(s)  =  - 

(  0)(  .7,  .25] (  5) 


.  222 

YW ( s)  =  - 

[  -7,  .5] 


Pilot  models  and  low  order  approximations 

202  (  0)  (  . 0299)  [  .7,  -  25 ] [  .71,  .296][  .7,  .5](  1.44)(  3.24)(  5) 

(  5.17)  (  10.1)  (  10.2)[-. 866,  17.3] 

VP(s)  =  — - - - - - - . — . — 

(  0) (  . 04 ) (  .7,  .25] [  .7,  . 5] [  .7,  .5) (  .9)(  3.37)(  5) (  10.1) 

(  10.1)(  .259,  12.2] [  .864,  24.6] 


2  02  (  . 0299 )  [  .71,  .296](  1.44)(  3.24)(  5)  (  5.17)(  10.2) 

[-.866,  17.3] 

YPl(s)  =  - 

(  . 04 )  [  .7,  . 5 ]  (  . 9 )  (  3.37)  (  5)  (  10.1)  [  .259,  12.2][  .864,  24.6] 


-  -  677  (  . 0299 )  [  .71,  .  296  )  (  1.44)  (  3.24)  (  5)  (  5.17)(-9.37) 

YP2 ( s )  =  - 

(  .04)  [  .7,  .5]  (  .9)  (  3.37)  (  5)  (  9.37) 


-2 . 61  (  1.44)  (  5.17)  (-9.37) 

YP3 ( s )  =  - 

(  -9) (  9.37) 


Figure  9:  Summary  of  Dander  Single  Axis  Example 
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coo.  O  a  C  *  3  (Q 


Figure  10:  Frequency  Responses  for  Dander  Single  Axis  Example 
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K  <1  _  h  nt  y  J 

1  .75  "l  .5  500  10"fi  27 


and  the  full  attention  case  for  the  .(//  axis  (axis  #3)  is  summarized: 


Vc(s) 


10(0.1  i  .53 

(3.0)[.5.  .5]  [.7.. 5] 


ly  t\  >i  o]  g  J 

1  .07  .14  .si  .0013  .05 


Each  of  the  single  axis  OC'M  problems  are  solved  for  several  different 
attentional  fractions.  How  to  do  this  using  Program  C  is  now  explained.  As 
seen  in  Table  2.  the  parameter  /  is  stored  in  the  1x6  element  of  the  P\2 
matrix.  For  each  different  value  of  /  change  the  P12  matrix  and  repeat  the 
OCMKBF  macro.  This  will  start  t  he  K  BF  iterations  at  the  last  values  of  the 
noise  intensities.  The  iteration  summaries  for  the  attentional  fraction  cases 
/=  1/2.  1/3,  1/1.  and  1/5  are  listed  in  Figure  11.  The  commands  used  for 
this  attentional  fraction  survey  are  listed  below  (see  also  the  EX3  macro): 


CC>pl2(l ,6)=l/2  ft  GOCMKBF , . 1 
CC>pl2(l  ,6)=l/3  ft  (OOCMKBF,.  1 
CC>pl2(l  ,6)  =  l/4  ft  (OOCMKBF,.  1 
CC>pl2(  1 ,6)  =  l/4  ft  (OOCMKBF,.  1 

For  each  axis  create  a  table  of  optimal  costs  .)  versus  attentional  fractions 
/.  Normalize  the  costs  by  dividing  by  mean  square  error  of  the  driving  noise. 
It  has  been  empirically  determined  that  a  linear  relationship  exists  between 
./  and  1//  over  a  range  of  attentional  fractions,  which  are  for  this  example: 


JiK 

=  .017//,  +  .028 

for  /,  >  .1 

■h/crf2 

=  .668/ /  -  .212 

for  f2  >  .5 

■h/crl 

=  -018//;,  +  .01  1 

for  /3  >  .13 

The  total  optimal  cost  is  the  sum  of  ./, .  ./2,  and  ./ 3.  The  following 
attentional  fractions  minimize  the  total  cost: 

1  /  f\  -  1  +  +  \J  f>i/«  l 

1  /  fi  =  IT  \J <i\  I  'i  2  +  \J «-j  / 

1  /  /a  =  1  +  \Ja\  / "  i  +  3 
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Iteration  #  5  f  =  1/2 

Optimal  cost:  q,  g,  output,  input  rate,  total 

1 . OOOOE+OO  1 . 8000E-04  7.2827E-03  2.0355E-03  9.3182E-03 

Performance  :  E{y_lA2),  E(y_2A2),  E(u_aA2),  E{uA2),  E{(du/dt)A2) 

7 . 2827E-03  5.8567E-02  2.2780E-01  1.8749E-01  1.1308E+01 

Old  noise  intensities  (V_yl,  V_y2 ,  V_ua) :  6.9377E-04  4.6353E-03  7.0520E-03 

New  noise  intensities  (V_yl,  V_y2 ,  V_ua) :  7.0840E-04  4.7101E-03  7.1565E-03 

Noise  ratios  dB  (rho_yl,  rho_y2,  rho_ua) :  -20.0907  -20.0695  -20.0639 

Max  noise  ratio  difference  =  9.065247E-02  dB,  Threshold  =  .1  dB 

Iteration  #  9  f-1/3 

Optimal  cost:  q,  g,  output,  input  rate,  total 

1 . OOOOE+OO  1 . 8  000E-04  9.7001E-03  2.2943E-03  1.1994E-02 

Performance  :  E{y_lA2),  E{y_2A2),  E(u_aA2),  E(uA2),  E((du/dt)A2) 

9.7001E-03  6.7465E-02  2.5446E-01  2.0850E-01  1.2746E+01 

Old  noise  intensities  (V_yl,  V_y2,  V_ua) :  1.3042E-03  7.8898E-03  7.8965E-03 

New  noise  intensities  (V_yl,  V_y2 ,  V_ua) :  1.3288E-03  7.9950E-03  7.9940E-03 

Noise  ratios  dB  (rho_yl,  rho_y2 ,  rho_ua) :  -20.0813  -20.0575  -20.0533 

Max  noise  ratio  difference  =  3.126068E-02  dB,  Threshold  =  . 1  dB 

Iteration  #  13  f-l/d 

Optimal  cost:  q,  g,  output,  input  rate,  total 

1. OOOOE+OO  1 . 8000E-04  1.2166E-02  2.5227E-03  1.4689E-02 

Performance  :  E(y  1A2),  E{y_2A2},  E{u_aA2),  E(uA2),  E{(du/dt)A2} 

1 . 2166E-02  7 . 5342E-02  2.7794E-01  2.2698E-01  1.4015E+01 

Old  noise  intensities  (V_yl,  V_y2 ,  V_ua) :  2.0933E-03  1.1619E-02  8.6402E-03 

New  noise  intensities  (V_yl,  V_y2,  V_ua) :  2. 1288E-03  1.1752E-02  8.7318E-03 

Noise  ratios  dB  (rho_yl,  rho_y2,  rho_ua) :  -20.0731  -20.0491  -20.0458 

Max  noise  ratio  difference  =  7.305717E-02  dB,  Threshold  =  . 1  dB 

Iteration  #  17  f=  1/5 

otimal  cost:  q,  g,  output,  input  rate,  total 

1. 0000E+00  1 . 8000E-04  1.4690E-02  2.7294E-03  1.7419E-02 

Performance  :  E(y_lA2),  E{y_2'2),  E{u_aA2),  E(uA2),  E((du/dt)A2) 

1 . 4  690E-02  8.2502E-02  2.9919E-01  2.4369E-01  1.5164E+01 

Old  noise  intensities  (V_yl,  V_y2,  V_ua) :  3.0646E-03  1.5767E-02  9.3124E-03 

New  noise  intensities  (V_yl,  V_y2 ,  V_ua) :  3.1121E-03  1.5924E-02  9.3992E-03 

Noise  ratios  dB  (rho_yl,  rho_y2,  rho_ua) :  -20.0668  -20.0430  -20.0403 

Max  noise  ratio  difference  =  .0668335  dB,  Threshold  =  .1  dB 


I  igure  11:  Attentional  Fraction  Survey  for  Dander  Multi  Axis  Example 


The  results  for  this  example  are: 


(.12  .76  .12) 
1.03 


(  /i  fi  fz)  - 
J  = 

Considerably  more  cost  is  required  to  control  axis  #2,  as  might  be  expected 
from  the  unstable  dynamics  in  that  axis,  and  as  a  result  considerably  more 
attention  is  given  to  axis  #2. 

The  optimal  costs  just  computed  can  be  used  to  make  the  following 
prediction  of  pilot  opional  ratings  (PORs): 

HOR  =  5.5  +  3.7 log, o  f  ~y~j) 

J  is  the  optimal  cost  as  computed  using  an  error  weight  q=\,  <J2C  is  the  mean 
error  of  the  driving  noise,  is  the  square  of  the  input  noise  bandwidth,  and 
the  numbers  5.5  and  3.7  are  experimentally  determined  using  the  Dander 
data  and  the  new  STI  data  in  this  report.  A  ±1  spread  in  the  in  the  PORs 
is  tolerable.  The  results  for  this  example  are: 


POR  POR 


axis 


J 


„2 


(predicted)  (experimental) 


1 

.015 

2.7 

2,5-3 

2 

.156 

6,1 

5  0.5 

3 

.062 

3.2 

3  3.5 

all  3 

1 .03 

7.7 

7-C 

The  predictions  in  this  case  are  very  good.  All  of  the  9  different  cases 
of  the  3  axis  Dander  experiments  have  been  similarly  analyzed  using  the 
OCM  implemented  in  Program  CC.  Predicted  PORs  were  within  ±1  of 
the  experimental  results  in  7  out  of  9  cases.  The  POR  predictions  for  the 
Dander  data  agrees  case-by  case  with  the  analysis  reported  by  McRuer  and 
Schmidt.  [16]. 


A  How  to  Use  Program  CC,  Version  4 

Program  CC  is  a  command  driven  computer-aided-control-system-dosign 
package.  It  is  written  in  compiled  BASIC  and  operates  under  the  DOS 
operating  system  on  the  IBM- PC  and  compatible  personal  computers.  Pro¬ 
gram  CC  works  with  linear  systems,  either  analog  or  digital,  which  are 
modeled  with  either  transfer  functions  or  state  space  equations.  A  large 
number  of  classical,  sampled-data,  and  state  space  algorithms  are  imple¬ 
mented,  including  an  extensive  amounts  of  interactive  graphics.  The  user 
environment  in  Program  CC  is  robust,  friendly,  and  very  powerful. 

This  appendix  gives  an  overview  of  Program  CC  which  is  tailored  to  users 
of  the  human  optimal  control  model.  It  is  recommended  that  you  operate 
the  program  while  reading  this  introduction.  Execute  the  commands,  or  at 
least  some  of  them,  when  they  are  introduced. 

A.l  The  Basics 

There  are  about  300  commands  arranged  in  a  hierarchy.  Use  the  manual  or 
the  on-line  help  (command  HELP)  for  a  list  of  the  commands  and  parame¬ 
ters.  The  hierarchical  levels  of  interest  here  are: 

CC  transfer  function  commands 

STATE  state  space  commands 

DATA  data  file  commands 

MACRO  creating  and  editing  macros 

Commands  are  entered  in  response  to  a  prompt,  for  example: 

COHELP 

STATE>HELP 

Commands  and  parameters  are  entered  in  either  upper  or  lower  case 
letters.  Blanks  are  ignored,  and  only  the  minimum  number  of  letters  to 
make  a  command  non  ambiguous  is  needed.  Separate  the  parameters  by 
commas,  include  them  on  the  command  line  (expert  mode)  or  let  yourself  be 
prompted  (novice  mode).  Include  more  than  one  command  on  a  command 
line  by  separating  them  with  &. 

Branch  to  different  levels  of  the  command  hierarchy  by  typing  the  name 
of  the  command  level.  Those  users  familiar  with  Version  3  should  note 
that  the  BUILD  and  MR  command  levels  have  been  compressed  into  the 
CC  level.  Return  to  the  CC  level  by  typing  either  QUIT  of  CC.  Leave  the 
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program  from  the  CC  level  by  typing  either  QUIT  (prompts  with  Are  you 
sure?)  or  EXIT  (no  prompt). 

BASIC  execution  errors  are  trapped  by  the  program,  resulting  in  an  error 
message  and  a  return  to  the  CC  command  level.  The  following  commands 
are  built  into  the  program: 

Ctrl-NumLock  temporary  halt 

FI  halt  and  return  to  the  CC  command  level 

|  i  recall  previous  commands 

Use  the  Fl  key  to  abort  commands  and  macros.  If  waiting  for  a  prompt 
then  follow  Fl  with  a  carriage  return.  The  program  does  not  respond  to 
Fl  if  it  is  pressed  during  an  overlay  swap,  so  sometimes  it  has  to  be  pressed 
several  times. 

All  transfer  function  and  state  space  calculations  are  computed  using 
double  precision  arithmetic,  which  has  approximately  18  decimal  places  of 
accuracy.  The  data  files  are  stored  in  single  precision,  mainly  to  save  disk 
space  when  extreme  accuracy  is  not  needed. 

Version  -1  of  Program  CC  is  compatible  with  the  EGA  and  VGA  mon¬ 
itors.  but  still  retains  downward  compatibility  with  the  CGA.  The  high¬ 
est  screen  mode  is  automatically  determined  when  the  program  starts. 
Monochrome  monitors  are  not  recommended  because  they  do  not  allow 
graphics,  though  state  space  operations  (in  fact  to  entire  OCM  solution) 
are  still  possible.  Switch  screen  modes  with  the  SCREEN  command.  Four 
color  plots  can  be  created  on  the  EGA  and  VGA.  Switch  colors  with  the 
command  COLOR. 

Program  CC  is  made  of  many  separate  modules,  which  are  chained  as 
needed  to  execute  commands.  The  modules  and  the  data  can  be  located  on 
more  than  one  disk,  as  determined  by  pointers  established  in  the  SETUP 
command.  Placing  the  data  and  the  most  used  modules  in  a  RAM  disk 
significantly  speeds  up  execution.  User  defined  commands  can  be  included 
in  the  list  of  '('  commands,  again  as  established  in  the  SETUP  command. 
When  the  user  defined  command  is  called  then  the  specified  program  is 
chained.  All  of  the  information  established  by  the  SETUP  command  is 
stored  in  the  SSSETUP  file,  and  is  therefore  available  each  time  the  program 
starts.  See  Section  3.2.1  for  information  on  the  best  setup  for  the  OCM 
macros. 

There  are  several  more  miscellaneous  commands  of  general  interest: 
FILES  lists  directory  files,  NAME  renames  them,  and  KILL  kills  them. 


CLS  clears  the  screen.  TIME  and  DATE  provide  their  namesakes.  CLOCK 
draws  an  analog  clock.  CALCULATOR  is  an  HP-like  calculator.  ECHO 
echoes  a  message  and  PAUSE  pauses  until  a  key  is  pressed,  both  of  which 
are  useful  in  macros.  SHELL  shells  to  DOS,  and  EXIT  returns  to  CC. 
Use  DESQVIEVV  to  run  multiple  copies  of  CC,  or  CC  and  any  other  set 
of  programs,  and  switch  between  the  programs  with  simple  DESQVIEVV 
commands. 

A. 2  Data  Types 

There  are  four  data  types: 

Transfer  function 
State  space  quadruple 
Transfer  function  matrix 
Data  file 

The  names  can  be  six  letters  or  numbers,  starting  with  a  letter,  followed 
by  a  prefix  of  3  letters  or  numbers.  For  example  G,  P,  YP ,  7/100,  and 
YP.G.  The  convention  in  earlier  versions  of  Program  CC,  which  is  still 
largely  maintained  in  the  documentation,  is  to  use  the  names  G,  for  trans¬ 
fer  functions,  Pt  for  state  space  quadruples,  and  //,  for  transfer  function 
matrices  (the  subscript  ?  refers  to  any  positive  integer). 

ASCII  data  files  are  automatically  created  on  disk  for  each  data  element. 
The  files  are  easily  recognized  because  they  always  begin  with  $$.  The  files 
remain  on  the  disk  after  the  program  is  finished,  do  matter  how  ungracefully, 
and  are  therefore  available  when  the  program  is  restarted. 

A  transfer  function  is  a  ratio  of  polynomials.  The  polynomials  can  be 
factored  in  any  desired  way.  A  state  sjhicc  quadruple  is  a  packed  set  of 
four  real  matrices  which  represent  a  state  space  differential  equation.  Real 
matrices  are  a  subset  of  state  space  quadruples,  being  just  the  constant  term. 
A  transfer  function  niatrii  is  a  matrix  of  transfer  functions,  usually  obtained 
by  converting  from  a  multivariable  state  space  quadruple.  Finally,  a  data 
file,  is  an  indexed  set  of  real  or  complex  matrices  where  the  index  represents 
either  times  or  frequencies,  and  the  matrices  are  ordered  row  wise. 

The  systems  are  either  analog  or  digital,  depending  on  a  flag  set  with  the 
commands  ANALOG  and  DIGITAL.  The  program  starts  with  ANALOG  as 
the  default.  The  OCM  deals  only  with  analog  systems,  and  therefore  the 
sampled-data  features  are  not  emphasize  in  this  introduction. 
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A. 3  The  CC  Command  Level 


Transfer  functions  can  be  entered  either  by  coefficients,  by  shorthand  nota¬ 
tion,  or  symbolically.  For  example,  to  enter: 


G(s)  = 


10(s+  1) 


s[s2  +  2(.l)(10)s  +  (10)2] 

use  one  of  the  following  commands: 

COGENTER.g,  2,0,10,1,1,1,  2,1,1,0,2,1,2,100 
CC>SENTER,g ,  2,0,10,1,1,  2 , 1 ,0 ,2 , . 1 , 10 
CC>g=10*(s+l)/s/(s*2  ♦  2*. 1*10*8  +  10*2) 

Let  yourself  be  prompted  for  the  coefficients  until  you  understand  their  or¬ 
der.  Most  users  prefer  to  use  the  symbolic  form.  Several  commands  are  avail¬ 
able  for  particular  types  of  transfer  functions:  BUTTERWORTH,  CHEBY- 
SHEV,  BESSEL.  PADE.  LEADLAG.  INTEGRATOR.  NOTCH,  and  soon. 

Individual  coefficients  can  be  changed  using  the  equation  interpreter. 
Its  best  to  precede  a  change  with  a  display  of  the  transfer  function.  In  the 
following  example,  the  denominator,  2nd  polynomial,  Oth  order  coefficient 
is  changed  to  200: 

CODISPLAY,  g 
CC>g(d,2,0)=200 

Transfer  functions  can  be  displayed  in  many  different  forms,  as  sug¬ 
gested  by  the  command  names:  DISPLAY,  SHORTHAND,  SINGLE,  and 
1' NTTARY.  PZF  (pole  zero  form).  TCF  (time  constant  form).  PFE  ( part ial 
fraction  expansion),  and  ILT  (inverse  Laplace  transform).  In  the  C’C  com¬ 
mand  level  entering  just  the  name  is  equivau.it  to  DISPLAY.6',.  Try.  foi 
example: 

COG  ft  SINGLE  ft  PZF  ft  SH0  ft  PFE  ft  ILT 


The  transfer  function  coefficients  are  stored  in  as  shown  by  the  DISPLAY 
command.  Change  the  polynomial  factors  with  the  commands  CHPZF, 
ClITCF.  CIISINGLE,  and  CH UNITARY.  For  example: 

COG  ft  CHSINGLE,G,G1  ft  G1 

Transfer  functions  can  be  built  up  from  others  using  the  equation  in¬ 
terpreter.  This  powerful  facility  is  used  by  Program  CC  in  place  of  block 
diagrams  to  compute  loop  and  closed  loop  transfer  functions.  For  example: 
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CC>G1=G*G1 

CC>G3= ( G+Gl ) / ( 1+G2* (G+Gl)  ) 

CC>G10=(s+l)*G/s 

The  following  example  is  a  simple  way  to  find  the  closed  loop  poles  with 
unity  feedback: 

CC>G2=G/ (1+G)  4  PZF.G2 

The  same  can  be  done  with  the  command  STABILITY.  The  range  of  stable 
gain  can  be  found  with  the  command  ROUTII. 

The  equation  interpreter  cancels  polynomial  factors  if  the  coefficients 
are  linear  multiples,  but  factors  are  not  further  broken  down  to  look  for 
cancellations.  For  example,  (2s  +  2 )/(s  -t-  1)  is  cancelled  but  not  (s2  +  2s  + 
l)/(s  +  1).  The  NEAR  command  converts  to  first  and  second  order  factors 
and  cancels  with  e  tolerance. 

The  displays  default  to  single  precision.  Change  this  to  double  or  mini 
(3  decimal  places)  precision  using  the  FORMAT  command.  Line  printer 
output  of  transfer  functions  is  produced  by  preceding  the  display  style  by 
LP.  For  example:  LPDISPLAY,  LPSIIORTHAND,  LPPZF,  and  LPTCF. 
Messages  can  be  included  on  the  output  by  surrounding  the  message  with 
quotes,  for  example: 

COLPSHQ.Gl, "Closed  loop  system" 

Transfer  functions  can  be  stored  and  recalled  under  arbitrary  file  names 
with  the  commands  STORE  and  RECALL,  which  is  useful  for  long  term 
saving  of  special  files.  Large  numbers  of  transfer  functions  are  more  eco¬ 
nomically  saved  by  using  DOS  to  copy  them  to  a  sub-directory. 

An  e“ai  delay  can  be  included  with  the  transfer  function  using  the  DE¬ 
LAY  command.  All  transfer  functions  are  multiplied  by  the  same  delay,  and 
transfer  functions  with  delays  cannot  be  algebraically  combined.  It  is  useful, 
nevertheless,  to  use  exact  delays  for  frequency  response  calculations. 

Still  more  transfer  function  commands  are  available.  ADJOINT  com¬ 
putes  either  G,[-s)  or  G,(l/z).  PARTIAL  extracts  terms  from  a  partial 
fraction  expansion,  and  SPECTRAL  extracts  poles  in  a  half  plane.  MEAN- 
SQUARE  computes  the  mean  square  error.  NEAR  cancels  poles  and  zeros 
according  to  either  an  absolute  or  relative  criterion.  LFAPPROX  replaces 
high  frequency  poles  and  zeros  with  a  Pade  delay  approximation,  and  HFAP- 
PROX  does  something  similar  with  low  frequency  poles  and  zeros.  Use  the 
equation  interpreter  to  extract  numerators  and  denominators  and  to  make 
substitutions  for  .s: 


CC>nd=G(n)/Gl(d) 

CC>gminus=g(-s) 

CC>gf-g(f) 

The  following  commands  require  more  explanation  than  is  appropriate 
here:  DIOPOLE  and  DIOLQG  are  pole-placement  and  LQG  algorithms 
solved  by  Diophantine  equations.  INNER,  OUTER,  BLASHKY,  WIENER, 
FILTER,  and  LQG  are  Wiener-llopf  commands. 


A. 4  The  STATE  Command  Level 

State  space  quadruples  are  entered  using  the  command  PENTER,  which 
prompts  for  the  number  of  states,  inputs,  and  outputs;  and  then  prompts 
for  the  respective  row-wise  elements  of  the  A,  B,  C,  and  D  matrices.  All  of 
this  can  be  placed  on  the  command  line,  for  example: 

STATE>PENTER , P , 2 , 1 , 1 ,  1,1, 1,1,  2,2,  3,3,  4 


creates  the  quadruple: 


/I 

1 

2 

1 

1 

2 

\  3 

3 

1 

which  represents  the  state  space  differential  equation: 


dt 


u 


Real  matrices  are  lust  the  D  part  of  the  quadruple.  The  following  command: 
STATE>PENTER ,P1 ,0,2,3,  1,2, 3, 4, 5, 6 


creates  the  real  matrix: 

1  2 
3  1 

5  6 

Real  matrices  can  be  entered  symbolically  using  the  state  space  equa¬ 
tion  interpreter.  Elements  in  a  row  are  separated  by  commas  and  rows  are 
separated  by  semicolons.  Enter  and  then  display  the  previous  P\  matrix  as 
follows: 
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STATE>P1=(1 ,2; 3, 4; 5,6)  ft  Pi 


The  equation  interpreter  does  not  allow  direct  entry  of  state  space  quadru¬ 
ples,  but  they  can  first  be  entered  as  a  real  matrix  and  then  partitioned  by 
changing  the  number  of  states: 

STATE>P=(1,1,2;  1,1,2;  3,3,4)  ft  P=CHST(P,2)  ft  P 

The  command  DISPLAY,/5,  (or  just  Px )  is  used  to  display  Px.  To  display 
only  part  of  the  quadruple  use  the  command  P(A).  To  display  the  elements 
in  double  precision  use  the  FORMAT  command.  To  display  just  the  di¬ 
mensions  use  WHAT, Pi.  To  eliminate  elements  with  absolute  values  less 
than  a  threshold  use  Pl  =  EPS( PI )-  To  obtain  a  line  printer  listing  use  the 
LPDISPLAY,P1  command,  which  can  contain  a  message,  e.g.: 

STATE>LPDIS ,P15 , "KBF  gains" 

Single  elements  of  a  matrix  or  a  quadruple  can  be  changed  by  referring 
to  their  position: 

STATE>P1(1,2)»10 

STATE>P(C , 1 ,2) =100 

This  is  the  easiest  type  of  augmentation.  More  generally,  a  real  matrix 
can  be  augmented  at  the  >,jtlx  position,  which  overwrites  a  block  of  the  old 
matrix  and  extends  the  boundaries  if  necessary. 

Real  matrices  (and  most  quadruples)  can  be  nested  and  built  up  from 
smaller  matrices.  For  example,  to  add  a  row  to  Pw: 

STATE>P10=(P10 ;  1,1) 

The  new  matrix  is  built  up  using  augmentation  with  zero-fill,  so  there  is 
never  a  problem  with  imsn,<uc!ieu  numoers  ol  rows  and  columns,  lo  save 
some  typing,  for  example,  the  3x3  identity  matrix  can  be  entered  in  the 
following  way: 

STATE>p3=(l ;  0,1;  0,0,1) 

Parts  of  a  real  matrix  or  a  state  space  quadruple  can  be  extracted,  as 
shown  in  the  examples: 
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STATE>P=P1(1 ,2)  single  element  of  P\ 

STATE>P=P1(1)  1st  row  of  Pi 

STATE>P=P1  (  ,2)  2nd  column  of  P, 

STATE>P=P1  ( 1 :3 ,4: 5)  rows  and  columns  of  P\ 

STATE>P=P1  (A)  P\[A)  from  the  Pi  quadruple 

STATE>P=P1(C,  ,4:5)  columns  4  and  5  of  PX{C) 

Several  special  types  of  matrices  constructed  with  the  following  func¬ 
tions:  IDEN,  DIAG,  and  RND.  Controllable  canonical  realizations  of  stan¬ 
dard  filters  are  created  by  the  commands  BUTTERWORTH,  CHEBYSHEV, 
BESSEL,  PADE.  LEADLAG,  INTEGRATOR,  and  NOTCH,  and  so  on. 
Transfer  functions  can  be  converted  to  state  space  quadruples  with  the  com¬ 
mands  CCF,  OCF,  and  DCF;  respectively  controllable,  observable,  and  di¬ 
agonal  canonical  forms.  Transfer  functions  included  in  state  space  equations 
are  automatically  converted  using  the  CCF  realization.  Try  the  following: 

CC>g=10+(s+l)/(s~3  +  2*s*2  +  3*s+4) 

COSTATE 

STATE>p=g 

STATE>OCF,g,pl 

The  number  sign,  #,  is  used  for  diagonal  augmentation.  It  is  used  in  the 
following  commands  to  create  2nd  order  Fade  approximations  of  diagonal 
delays: 


STATE>PADE,P10, .1,2  &  PADE, Pll , .2,2 
STATE>DELAY=(P10#P11) 

The  state  space  equation  interpreter  can  be  used  to  algebraically  combine 
state  space  quadruples.  Dimensions  are  checked  for  validity,  and  when  state 
space  quadruples  are  used,  the  appropriate  combinations  of  the  A ,  B,  C, 
and  D  matrices  are  automatically  computed.  The  valid  operations  are: 

4-  Addition 
*  Multiplication 
Exponentiation 
'  Transpose 

/  Right  division,  a/i3-n  *  ;3~] 

\  Left  division,  a\J  —  a'1  ♦  J 
|  Feedback,  o|J  =  a  *(  1  +  ,3  +  o  )-1 
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The  levels  of  precedence  are  as  follows,  with  operations  on  the  same  level 
computed  left  to  right  as  they  appear  in  the  equation: 

(  ) 

t 

*/\\ 

+  - 

An  mnemonic  for  left  and  right  multiplication  is  that  the  top  of  the  slash 
points  to  the  inverted  quadruple.  The  |  operation  for  feedback  results  in  a 
minimal  order  state  space  realization,  with  a  in  the  feedforward  and  (3  in  the 
feedback  paths.  Identity  matrices  in  the  equation  are  automatically  sized. 
Try,  for  example,  the  following: 

STATE>BUTTER,P , 1 ,4 
STATE>P1=P/(I+P) 

STATE>P2=P| I 
STATE>WHAT,P1  ft  WHAT.P2 

Both  F]  and  F2  are  closed  loop  systems  with  unity  feedback,  but  P\  is  a 
non-minimal  8th  order  realization,  and  P2  >s  the  much  preferred  minimal 
4th  order  realization. 

There  are  several  ways  to  invert  a  matrix  (a  quadruple  is  invertible  if 
the  D  term  is  invertible): 

STATE>P1=I/P 

STATE>P1=P\I 

STATE>P1=P*-1 

Gaussian  elimination  (LU  decomposition)  with  scaling  and  partial  pivoting 
is  used  to  compute  the  matrix  inversion,  which  takes  only  order  n3  op¬ 
erations.  the  same  as  matrix  multiplication.  A  warning  is  printed  if  the 
determinant  is  zero.  It  is  more  efficient  not  to  invert  an  entire  matrix,  with 
the  first  example  below  being  preferred  over  the  second: 

STATE>P2=P1\P 
STATE>P2=Pl\I  ft  P2=P2*P 

The  PACK  command  packs  four  real  matrices  of  compatible  dimensions 
into  a  state  space  quadruple,  and  the  UNPACK  command  does  the  converse. 
In  the  following  example  the  F  quadruple  is  unpacked  and  a  second  ,  of 
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outputs  is  c  reated  which  are  the  derivative  of  the  former.  This  is  only  valid 
if  P{D)= 0.  The  derivatives  are  created  using  the  identity: 

Cs(sl~  A)~'B  =  CB  +  CA(sI  -  A)-1  B 

STATE>UNPACK,P, P500 ,P501 ,P502 ,P503 
STATE>P502* (P502 ; P502*P500) 

STATE>P503= (P503 ; P502*P50 1 ) 

STATE>PACK,P500 ,P501 ,P502 ,P503 ,P 

A  complicated  system  can  be  built  up  using  the  state  space  equation 
interpreter.  Series  elements  are  multiplied;  parallel  elements  are  either  added 
or  augmented,  depending  on  the  input  and  output  connections;  and  feedback 
paths  are  closed  using  the  |  operation.  Use  the  FEEDBACK  command  to 
feedback  a  subset  of  outputs  to  a  subset  of  inputs.  In  the  following  example, 
define  P  to  have  4  inputs  and  4  outputs,  and  define  P\  to  be  SISO.  A  single 
feedback  loop  is  closed  around  P  from  the  2nd  output  to  the  3rd  input,  with 
Pi  in  the  feedback  path: 

STATE>P2=P1 I ( (0 ; 0 ; 1 ; 0) *P1* (0 , 1 , 0 , 0) ) 

Functions  are  used  in  equations.  All  of  the  usual  trigonometric  and  expo¬ 
nential  functions  are  available,  and  when  applied  to  matrices  work  element- 
by-element.  Other  functions  which  are  available  include  DET,  TRACE, 
and  NORM.  Most  of  the  complicated  operations,  however,  are  computed 
using  commands.  The  NORM  command  computes  6  different  matrix  norms. 
PSEUDO  INVERSE  uses  the  singular  value  decomposition  (SVD)  to  com¬ 
pute  the  pseudo  inverse  of  a  rectangular  matrix.  SPACE  computes  orthonor¬ 
mal  basis  vectors  for  the  fundamental  subspaces.  PROJECTION  computes 
projection  matrices  onto  the  fundamental  subspaces. 

Several  different  matrix  decompositions  can  be  computed:  EIGEN¬ 
VALUE.  SCHUR.  HESSENBERG.  and  SVD.  Place  a  G  in  front  of  the 
first  three  commands  to  compute  generalized  versions.  The  following  exam¬ 
ple  computes  the  eigenvalue  decomposition  and  the  verifies  the  result.  The 
eigenvalues  are  stored  in  a  2  column  matrix,  which  is  converted  to  block 
diagonal  with  the  DIAG  function. 

STATE>EIG ,p ,d , x 
STATE>p-x*DIAG(d)/x 

Several  different  matrix  equations  can  be  solved:  LYAPUNOV,  RIC- 
C  All.  and  S't  L\  LSTOR.  Place  a  D  in  front  of  the  first  two  commands  to 
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compute  discrete  versions.  Several  different  methods  are  available  to  solve 
these  equations,  with  the  defaults  methods  being  recommended. 

The  SIMILARITY  command  computes  similarity  transformations,  ei¬ 
ther  with  respect  to  an  input  matrix,  or  with  respect  to  eigenvectors,  Schur 
vectors,  or  Hessenberg  vectors  of  the  system  matrix.  The  CONTROLLA¬ 
BILITY  and  OBSERVABILITY  commands,  which  are  identical,  each  check 
for  both  uncontrollable  and  unobservable  modes;  but  the  algorithm  is  not 
foolproof,  and  if  in  doubt  compute  a  similarity  transformation  with  respect 
to  the  eigenvectors,  and  check  for  zero  rows  of  B  for  uncontrollable  modes 
and  zero  columns  of  C  for  unobservable  modes. 

Conversion  from  state  space  to  transfer  function  matrices  is  computed 
using  commands  FADEEVA  and  CIEP.  The  former  uses  the  Fadeeva  algo¬ 
rithm,  which  is  very  fast,  but  computes  polynomial  coefficients  and  is  there¬ 
fore  unreliable  for  large  systems  (>  6th  order  depending  on  the  dynamics). 
The  Fadeeva  algorithm  is  used  in  the  following  example: 

STATE>FADEEVA,P,H 

STATE>EXTRACT , H , G , 1 , 1  k  EXTRACT, H, G1 , 1 ,2 

The  transfer  function  matrix  H(s)  is  created,  and  then  the  EXTRACT 
command  is  used  to  extract  elements  of  Il{s)  into  transfer  functions.  The 
command  EXTRACT, II, G, ALL  extracts  all  of  the  elements  of  H($)  row¬ 
wise  into  G(s).  G\ (s),  and  so  on.  Display  transfer  function  matrices  with 
the  same  commands  used  for  transfer  functions.  The  denominator,  which  is 
common  to  all  of  the  elements,  is  only  displayed  once. 

The  GEP  command  uses  the  generalized  eigenvalue  problem  (hence  the 
name)  to  compute  the  poles  and  zeros,  which  is  slower  but  numerically  more 
reliable.  One  problem  with  the  GEP  command  is  that  it  has  trouble  with 
zeros  at  infinity  (1/s2  has  2  zeros  at  infinity),  and  for  reasons  known  only 
to  a  few  tends  to  put  them  at  10Hl/,m.  where  m  is  the  order  of  the  infinite 
zeros.  The  GEP  commands  allows  the  user  to  set  a  magnitude  threshold, 
above  which  the  zeros  are  considered  infinite.  If  the  state  space  system  is 
SISO  then  the  conversion  can  be  directly  placed  into  a  transfer  function,  as 
done  in  the  following  example,  which  has  set  a  threshold  of  1012  for  infinite 
zeros: 


STATE>GEP ,P1 ,G1 , lel2 

The  POLE  PLACEMENT  command  uses  '.he  Ackermann  formula  for 
SISO  state  space  pole  placement  designs.  In  the  following  example,  the  pole 
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placement  method  is  used  to  compute  full  state  feedback  gains  (prompts  not 
shown),  and  then  the  closed  loop  system  is  computed  in  two  different  ways: 

STATE>POLE  PLACEMENT.P.Pl,  ... 

STATE>P2=P  ft  P2(A)=P2(A)-P2(b) *P1 
STATE>FEEDB ACK , 4 , P , P 1 , P2 

Optimal  control  designs  are  solved  using  the  commands  LQR,  KBF,  and 
LQG.  Digital  versions  are  solved  with  the  commands  DLQR.  DKBF,  and 
DLQG.  If  the  intermediate  Riccati  solution  is  desired,  it  must  be  obtained 
with  a  separate  call  to  the  RICCATI  command.  Several  different  options  are 
available  for  entering  the  data,  and  several  choices  of  computation  schemes 
are  available.  Its  too  much  to  explain  here,  and  one  example  will  have  to 
suffice.  The  1st  option  of  the  LQR  command  is  used  in  the  following: 

STATE>LQR, 1 ,P ,Q ,R ,F , le-4 

where  the  parameters  are  respectively  the  system,  state  weight,  input  weight, 
and  LQR  gains.  The  input  weight  is  multiplied  by  10-4.  The  default  com¬ 
putational  method  using  Schur  vectors  is  used  to  solve  the  Riccati  equation. 
In  the  following  example.  G  is  the  Riccati  solution,  and  F  is  the  LQR  gain: 

STATE>F=P(A)  ft  G=P(B) / ( le-4*R) *P(B) 

STATE>RICCATI .F.Q.G.G 
STATE>F=(le-4*R)\P(B) ’ *G 

While  these  several  commands  are  not  too  complicated,  if  this  sequence  is 
to  be  repeated  t lion  a  macro  should  be  constructed,  as  shown  later. 

A. 5  GRAPHICS 

Graphics  is  one  of  the  best  features  of  Program  CC.  Frequency  and  time 
domain  plots  can  be  obtained  starting  from  transfer  functions,  state  space 
quadruples,  or  data  files. 

The  plots  can  be  interactively  changed,  a  valuable  user-friendly  feature. 
The  interaction  is  accomplished  using  ]>lot  option  blocks.  Figure  12  contains 
a  Bode  plot  of  G(s)  =  (s-f  .2)e~  3,/(.s+  10).  The  plot  option  block  is  located 
underneath  the  plot.  Press  ?  (or  /)  to  list  the  full  names  of  the  options,  as 
shown  in  the  bottom  of  Figure  12.  and  as  listed  below: 

A  Overplot  with  an  additional  line 
B  (  rosshair  cursor 
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C  Transfer  function  cursor,  includes  robustness  calculations 
D  Change  foreground  and  background  options 
E  Change  axis  limits,  includes  zoom 
F  Create  data  file 

G  Use  equation  interpreter  to  augment  transfer  function 
H  Hardcopy 

I  Automatically  fill  in  with  more  points 
J  Change  plot  colors 
K  Include  date  and  time  on  title 
L  Label  the  plot 

M  Manually  fill  in  with  more  points 
P  Replot 
Q  Quit 

S  Plot  magnitude  asymptote 
T  Toggle  thick  and  thin  lines 
W  Change  title 
X  Pause,  for  use  with  macros 
V  Change  how  phase  is  calculated 
Z  Center  plot 
?  Help 

The  following  plots  are  available  for  transfer  functions: 

BODE 

NICHOLS 

XYQITST 

LOG  XYQITST  (uses  log  axes) 

ROOT  LOCUS  (uses  gain  stepping  algorithm) 

FASTRL  (uses  curve  tracing  algorithm) 

SIGGY  (Siggy  and  Bode  root  loci  plots) 

TIME  (inverse  Laplace  transform) 

DTIME  (inverse  c-transform) 

SIMULATION  (simulation  of  Laplace  transform) 

DSIMULATION  (simulation  of  r  transform) 

The  selection  and  operation  of  the  plot  options  change  depending  on 
tlie  type  of  plot.  It  suffices  here  to  use  the  above  BODE  example.  The 
axis  limits  for  each  plot  can  either  be  manually  entered  or  automatically 
determined. 


Options)  ft  BCDEFGH  I  J  K  LHP  l)S  IUXVZ  ?=helP 


A-addline  E=chg  Units  I  =f ill  tt=nore  pts  T=thick  lines  Z=center 

B-crosshair  F^data  file  J=colors  P=replot  U^title  ?=belp 

C=tf  cursor  G=equation  K=clock  Q=quit  X-pause  tf=G 

D=chg  options  H^hardcopy  L= label  S^asynptote  ¥=chg  phase 


l  iguro  12:  Program  ('('  Plot  Including  the-  Plot  Option  Blocks 


The  frequency  plots  work  from  the  FREQ  data  file,  created  by  the  FRE¬ 
QUENCY  command.  The  FREQ  file  saves  time  when  switching  between 
different  frequency  files,  it  allows  non-uniform  spacing  of  points,  and  it  saves 
disk  space  because  the  file  is  overwritten  each  time  the  FREQUENCY  com¬ 
mand  is  called.  To  creat"  the  Bode  plot  in  Figure  12: 

CC>G=(s+.2)/(s+J0)  &  DELAY,. 3  t  SHO.G 

COFREQ.G,  .01,100,100 

CC>B0DE, 3, A 

As  usual,  novices  enter  just  the  command  names  and  let  themselves  be 
prompted.  Several  different  plot  options  were  used  to  obtain  the  finished 
plot. 

The  TIME  and  DTIME  commands  compute  functions  of  time  which 
are  then  plotted.  Options  exist  for  open  loop  impulse  and  step  responses, 
closed  loop  impulse  and  step  responses,  and  non-causal  step  responses  (e.g. 
for  autocorrelations).  Ramp  and  sinusoidal  responses  can  be  computed  by 
first  augmenting  the  transfer  function  and  then  plotting  an  impulse  response. 
The  following  example  plots  a  closed  loop  step  response  and  then  a  closed 
loop  ramp  response: 

COTIME, G,  1 , AUTO  (l=closed  loop  step) 

q 

CC>G1=G/(1+G)/S“2 

COTIME ,G1 ,4 , AUTO  (4=open  loop  impulse) 

The  SIMULATION  and  DSIMULATION  commands  compute  simula¬ 
tions  of  transfer  functions.  The  same  open  and  closed  loop  impulse  and  step 
response  options  are  available,  but  not  the  non  causal  option  Data  files 
can  be  used  as  input,  representing  either  an  output  from  another  system  or 
an  arbitrary  input  sequence  creat- *1  by  the  INPUT  command.  The  SIM¬ 
ULATION  command  combines  a  call  to  the  CONVERT  command  (zero 
order-hold  equivalence)  with  a  call  to  the  DSIMULATION  command.  As 
an  alternative: 

COCONVERT, G,G1, 3,  .1  (Bilinear  with  7  =.l) 

DIODSIM  ,G1 ,2  ,  AUTO  (2=open  loop  step) 

When  in  the  digital  model  the  C'C  command  level  prompt  changes  to  DIG. 
Switch  back  to  the  analog  mode  with  the  ANALOG  command. 

Time  and  frequency  plots  can  also  be  obtained  from  state  space  quadru¬ 
ples.  The  STATE  commands  SIMULATION  and  DSIMULATION  are  used 
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to  create  time  domain  data  files  (unlike  the  CC  commands,  then  do  not 
directly  result  in  a  plot).  The  simulations  can  be  multivariable,  and  several 
different  types  of  inputs  are  available,  including  arbitrary  data  files.  The 
SIMULATION  command  computes  a  matrix  exponentional  to  discretize  the 
system  and  then  performs  a  digital  simulation.  The  name  of  the  time  se¬ 
quence  file  which  is  produced  is  Pt.Y . 

The  FREQUENCY  command  in  the  STATE  command  level  is  used  to 
compute  a  frequency  data  file  for  a  multivariable  state  space  quadruple.  The 
name  of  the  output  frequency  file  is  the  system  name  with  the  file  prefix  .G. 

Data  files  are  plotted  with  the  PLOT  command.  Data  files  can  be  cre¬ 
ated  from  transfer  functions  using  plot  options,  by  using  the  DENTER  or 
INPUT  command,  or  as  just  described  by  several  STATE  commands.  More 
parameters  are  required  than  for  other  types  of  plots.  Data  files  are  indexed 
matrices  in  general,  therefore  choices  must  be  made  for  which  rows  and 
columns  to  plot  (ALL  for  everything).  Choices  must  also  be  made  for  what 
to  plot  on  the  horizontal,  the  left  vertical,  and  optionally  the  right  vertical 
axes.  The  axis  choices  are  T(time).  \V( frequency  rad/sec),  F(frcquency  Hz), 
R(real),  I( imaginary ),  M(magnitude),  and  P(phase).  Each  can  be  prefixed 
with  L(logio)  or  D(dB),  with  no  prefix  defaulting  to  linear.  For  example: 

TR  (time  plot  for  simulation  data) 

LWLM  (magnitude  Bode  plot  with  login  scaling) 

LNYLMP  (same  including  phase) 

IAVDM,  LWDMP,  LFDM.  WM.  FM.  WDM  (Bode  variations) 

R I  ( Nyquist ) 

PI.M.'PDM  (Nichols) 

In  the  following  example  a  state  space  simulation  is  computed  and  then 
plotted.  The  parameters  which  set  up  the  simulation  are  not  shown  here: 

STATE>SIMULATION ,P1 ,  ... 

STATE>PLOT , P 1 ,TR , AUTO 

And  a  multivariable  frequency  file  is  plotted: 

STATE>FREQ,P,  .01,100,60 
STATE>PLOT,  P.G,  LWLM,  AUTO 

I  he  FREQ  file  created  for  transfer  functions  can  be  directly  plotted  by  the 
PLOT  command,  though  the  plot  options  are  not  as  convenient: 

COFREQ.G,  .01,100,100 
COPLOT ,  FREQ ,  LWDMP  ,  AUTO 


A. 6  The  DATA  Command  Level 


Data  files  are  either  time  sequences  or  frequency  data  files  as  created  by 
Program  CC  commands.  The  file  structure  is  very  simple,  however,  and 
data  from  any  other  source  can  be  used.  The  file  structure  is: 

trows,  icolumns,  type  (0=real,  l=complex) 
index 

real  or  complex  matrix  entered  rowwise 
index 

real  or  complex  matrix  entered  row-wise 


The  numbers  are  stored  in  ASCII  format,  and  are  separated  by  blanks,  com¬ 
mas,  semi-colons,  or  carriage  returns.  Two  carriage  ret  urns  in  a  row  is  inter¬ 
preted  as  a  zero.  The  data  file  is  stored  on  disk  with  the  name  $$name.ext . 
and  referred  to  in  the  DATA  command  level  sans  $’s  as  name. ext. 

The  DATA  command  level  is  used  to  change  and  algebraically  combine 
data  files.  The  data  files  are  plotted  using  the  PLOT  command,  which  can 
be  called  from  any  command  level. 

An  equation  interpreter  can  be  used  in  the  DATA  command  level  to 
algebraically  combine  data  files,  for  example: 

DATA>t=pl -g*p -g 
DATA>t=t/(i+t) 

It  is  of  course  best  not  to  let  intermediate  files  proliferate. 

Among  the  most  sophisticated  commands  are  those  for  computing  eigen¬ 
values,  singular  values,  and  structured  singular  values  of  frequency  files: 
respectively  EIGENVALUE,  SVD.  and  EMU.  For  example,  to  obtain  a  sin¬ 
gular  value  plot: 

DATA>SVD ,  P.G,  P.S 

DATA>PL0T ,  P.S,  LWDM ,  ALL,  AUTO 

A. 7  Making  Macros 

There  is  a  great  reluctance  to  using  macros.  Somehow  they  must  oe  terribly 
difficult.  They  are  not  —  but  they  do  require  some  learning.  The  author 
freely  admits  reluctance  to  learn  macros  from  other  programs,  which  for  the 
most  part  are  no  more  difficult,  but  always  just  a  little  bit  different,  than 
those  in  Program  CC. 
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Follow  step  by  step  the  creation  of  the  STBODE  macro  for  creating  a 
standardized  Bode  plot: 

COMACRO 

MACROADD 

bode ,3,0, auto ,2,-60,60,6, auto 

MACROSTORE,  STBODE 

MACROQUIT 

COCSTBODE 

Now  create  a  different  version  of  the  STBODE  macro  which  first  com¬ 
putes  the  frequency  file  of  an  arbitrary  transfer  function: 

COMACRO 

MACR0>ADD 

freq.&l, .01,100,60 

bode ,3,0, auto ,2,-60,60,6, auto 

MACR0>ST0RE, STBODE 
MACR0>QUIT 
CC>«STB0DE, GS0 

The  DAMP  macro  creates  a  transfer  function  with  an  arbitrary  damping 
ratio,  after  first  echoing  a  message: 

COMACRO 

MACR0>ADD 

echo,  &1  =  transfer  function  with  &2  damping  ratio 
&l=10*(s+l)/s/(s"2+2*&2*10*s+100) 

MACR0>ST0R£,DAMP 

macro>;,;it 

COQDAMP  ,G10 , .  1 

I  lie  STATE  commands  RICCATI  and  I.QR  are  combined  in  the  fol¬ 
lowing  macro.  The  parameters  are  echoed,  and  then  the  PAUSE  command 
gives  the  user  to  abort  with  the  Fl  function  key. 

COMACRO 

MACROADD 

echo.  The  LQR1  Macro 


echo, 
echo , 
echo , 
echo , 
echo , 
echo , 
echo , 
pause 
state 


ftl  =  system 

62  =  state  weight  Q 

63  =  input  weight  R 

£4  =  input  parameter  rho 

65  =  Riccati  solution  P 

66  =  LQR  gains  F 

Hit  FI  to  abort,  any  other  key  to  continue 


&6=&1 (a)  &  &5=&1 (b) /(&4*&3)*&l(b) 
riccati ,&6,&2 ,&5 ,&5 
&6=(&4*&3)\ftl(b) ’  *&5 
echo.  End  of  LQR1  macro 


MACROSTORE, LQR1 
MACRO>QUIT 

COQLQR1  ,P,P1  ,P2  ,  le-4  ,P3  ,P4 

Okay,  so  t hat  one  wasn't  so  easy. 

The  MACRO  command  level  is  a  mini-text  editor.  So  far  you’ve  seen 
the  ADD  and  STORE  commands.  There  are  several  more  commands  to 
do  things  like  LIST.  SEARCH.  DELETE,  and  REPLACE:  but  rather  than 
explain  it  here  just  use  the  HELP  command. 

Everything  that  gets  entered  goes  into  HISTORY  file  (disk  file  S'SHIS): 
which  can  be  recalled  in  the  MACRO  command  level,  edited,  and  played 
back  as  a  demonstration.  Or  if  you  want  to  create  a  standard  plot  but  cannot 
remember  the  parameters,  then  create  the  plot  once  using  prompts  and  build 
a  macro  which  includes  the  prompts.  Plot  labels  get  stored  along  with  the 
cursor  tin  meats,  which  show  up  in  the  HISTORY  file  as  weird  symbols. 
Here's  a  simple  example,  with  no  editing,  of  a  demonstration  macro: 

COMACRO 
MACR0>HIST0RY 
MACR0>ST0RE , DEMO 
MACRO>qUIT 
COQDEMO 

If  your  disk  space  is  limited  and  the  HISTORY  file  is  getting  too  long,  delete 
it  by  leaving  CC  and  returning,  or  more  elegantly  by  using  the  TRUNCATE 
command: 


71 


COMACRO 

MACRO >TRUNCATE,0 
MACROQUIT 

Okay,  that’s  enough.  Macros  aren’t  so  tough.  I  won’t  pressure  you  into 
learning  more. 
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B  Macro  Listings 

OCMYCl 


echo,  OCMYCl  Macro 

echo,  Purpose:  Create  (Y_c)_all  with  driving  noise  at  output 
echo,  (Y_c)_all=(l;s)*(Gi,Gj) 
echo , 

echo,  SISO  controlled  system  (Gi)  =  ft  1 
echo,  SISO  noise  filter  (Gj)  *  &2 
echo,  MIMO  (Y_c)_all  (Pi)  =  &3 
echo , 

echo,  Hit  function  key  FI  to  abort,  any  other  key  to  continue 

pause 

state 

&3=(ftl ,&2) 

&3 (c ,0 , l)=&3(c) *(&3(a) ,ft3(b)) 
quit 

echo,  Finished  OCMYCl 


OCMYC2 

echo,  0CMYC2  Macro 

echo,  Purpose:  Create  (Y_c)_all  with  driving  noise  at  input 
echo,  (Y_c)_all=(l;s)*Gi*(l,Gj) 
echo , 

echo,  SISO  controlled  system  (Gi)  =  &1 
echo,  SISO  noise  filter  (Gj)  =  &2 
echo,  MIMO  (Y.c).all  (Pi)  =  &3 
echo , 

echo,  Hit  function  key  FI  to  abort,  any  other  key  to  continue 

pause 

state 

&3=&1*(1,&2) 

&3(c,0,l)=&3(c)*(ft3(a) ,&3(b)) 
quit 

echo,  Finished  0CMYC2 


OCMYC3 


echo,  0CMYC3  Macro 

echo.  Purpose:  Create  (Y_c)_all  with  general  driving  noise 

echo,  (Y_c)_all=(l;s)*Gi*(Gj ,Gk) 

echo, 

echo,  SISO  controlled  system  before  noise  (Gi)  *  41 
echo,  SISO  controlled  system  after  noise  (Gj)  *  42 
echo,  SISO  noise  filter  (Gk)  =  43 
echo,  MIMO  (Y.c).all  (Pi)  =  44 
echo , 

echo,  Hit  function  key  FI  to  abort,  any  other  key  to  continue 

pause 

state 

44=41* (42 ,43) 

44(c,0, l)=44(c)*(44(a) .44(b)) 
quit 

echo,  Finished  0CMYC3 


OCMALL 

echo,  OCMALL  Macro 

echo,  Purpose:  Complete  human  OCM  problem,  using  defaults 
echo , 

echo,  (Y_c)_all  =  41 
echo , 

echo,  Replaces  OCMLQR,  OCMSETUP,  and  OCMKBF 
echo,  FI  to  abort,  any  other  key  to  continue 
pause 

Qocmlqr ,41 , 1 , le-4 , . 1 , .001 

Socmsetup , 1 , . 2 , -20, -20 , -20 , le-2 , le-2,1,0,0,1 
Gocmkbf , . 1 
Oocmpilot ,.2,2 
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echo,  Finished  OCMALL 


OCMLQR 

echo,  OCMLQR  Macro 

echo.  Purpose:  Adjust  LQR  weights  for  tau_N 
echo , 

echo,  Y_c  *  41 
echo,  q  =  42 
echo,  g  =  43 
echo,  tau_N  =  44 
echo,  thresh  =  45 
echo , 

echo.  Use  Function  Key  FI  to  abort,  any  other  key  to  continue 

pause 

state 

p=41 

p2=42 

p3=&3 

quit 

ocm.lqr ,44,45 

echo,  Finished  OCMLQR  Macro 

echo,  Next  use  OCMSETUP  and  OCMKBF  Macros 


OCMSETUr 


echo , 
echo , 
echo , 
echo , 
echo , 
echo , 
echo , 


OCMSETUP  Macro 

Purpose:  Setup  KBF  iterations  for  Human  OCM 


intensity  of  driving  noise:  V_v  =  41 
visual  delay:  tau  =  42 
y.l  noise  ratio  (dB) :  rho_(y_l)  =  43 
y_2  noise  ratio  (dB) :  rho_(y_2)  =  44 


i  > 


echo,  u_a  noise  ratio  (dB) :  rho.(u.a)  =  45 

echo,  initial  y_l  noise  intensity  (non-zero):  V_(y_l)  =  46 

echo,  initial  y_2  noise  intensity  (non-zero):  V_(y_2)  =  47 

echo,  initial  u_a  noise  intensity  (zero  okay):  V_(u_a)  =  48 

echo,  y_l  indifference  threshold:  T.l  =  49 

echo,  y_2  indifference  threshold:  T_2  =  410 

echo,  axis  fractional  attention:  f  =  ft 1 1 

echo , 

echo,  Execute  after  QCMLQR 

echo,  Hit  Function  Key  FI  to  abort,  any  other  key  to  continuo 

pause 

state 

echo , - Augmentation 

p5=p4( 1 , cdim(p4) ) 
p6=(p4(l ,1 :cdim(p4)-l)/p5,0) 
p7=chst((-p5,p5; 1 ,0) ,1) 

echo, - Setup  KBF  problem 

p8=p ( , 1) *p7 

p9=p8(a)  4  p9=chst(p9 ,rdim(p9) ) 
expon ,p9 ,p9 ,&2  &  p9=chst(p9,0)  &  cinalog 
pl0»p(b,,2)*ftl*p(b,,2)» 
pi  1=0 

pl2s(&3,&4,&5, 49,410,41 1,42) 

pl3=(&6 ,47 ,48) 

quit 

echo,  Finished  OCHSETUP  Macro 
echo,  Next  use  OCMKBF  Macro 


OCMKBF 

echo,  OCMKBF  Macro 

echo,  Purpose:  KBF/linear  predictor  iterations  for  OCM 
echo , 

echo,  noise  ratio  threshold  (dB)  =  41 
echo , 

echo,  Execute  after  OCMLQR  and  OCMSETUP 

7<i 


echo,  Hit  Function  Key  FI  to  abort,  any  other  key  to  continue 
pause 

ocm.kbf ,&1 

echo,  Finished  OCMKBF  Macro 


OCMPILOT 

echo,  OCMPILOT  Macro 

echo,  Purpose:  Create  SISO  state  space  models  of  Y_p  and  Y_cl 
echo, 

echo,  Visual  delay;  tau  =  &1 

echo,  Order  of  Pade  approximation  of  delay  =  &2 
echo , 

echo,  Execute  after  OCMLQR ,  OCMSETUP,  and  OCMKBF 

echo.  Hit  function  key  FI  to  abort,  any  other  key  to  continue 

pause 

state 

echo, - Augmentation  preliminaries 

pade  ,p500  ,&1  ,&2 

p501  =  (p8(a) -pl5*p8(c) , -p8(b) *p500(d) *p6  ,p8(b) *p500(c) ) 

p502=(-p9*pl5*p8(c) , (p8(a)-p8(b)*p6; -p500(b)*p6,p500(a))) 

p501=(p501 ;p502) 

p502=(pl5;p9*p!5;0*p500(b)) 

p503=(0*p6,-p500(d)*p6,p500(c)) 

pack,p501 ,p502 ,p503 , ,p500 

echo , - Pilot 

p24=p500 

p24(b)  =  (p24(b, , l)+p24(a)*p24(b, ,2) ;p24(c) *p24 (b  ,  ,2) ) 
p24=p7  *p24( , 1) 

echo, - Closed  loop  system 

p25=(p* (p7* (p500 ,1)#1)) |-(1,0;0,1;0;0) 
p25=p25 ( 1 ) 

kill , $$p500  &  kill ,$$p501  &  kill,$$p502  &  kill,$$p503 
quit 

echo,  End  of  OCMPILOT  Macro 

echo,  Continue  analysis  using  OCMG  and  0CMFREQ1 


OCMG 

echo,  OCMG  Macro 

echo.  Purpose:  Convert  Y_p  and  Y_c  from  state  space  to  tfs 
echo, 

echo,  Y_p :  G_i  =  41 
echo,  Y_c:  G_j  =  42 
echo , 

echo,  Execute  after  OCMPILOT 

echo,  Hit  function  key  FI  to  abort,  any  other  key  to  continue 

pause 

state 

p500=p(l , 1) 
f adeeva,p500,42 
kill  ,$$p500 
gep ,p24 ,41 
quit 

near ,41 ,41,1 , le-4 
near ,42 ,42 , 1 , le-4 
echo,  End  of  OCMG  Macro 

echo,  Continue  analysis  with  regular  CC  commands 


OCMFKEQ1 

echo,  OCMFREQl  Macro 

echo,  Purpose:  Create  freq  plots  of  Y_p,  Y_c,  Y_p*Y_c,  Phi 
echo,  Assumes  Y_p=P24  and  Y_cl=P25  computed  using  OCMPILOT 
echo , 

echo,  low  freq  =  41 

echo,  high  freq  =  42 

echo,  #  points  (log  spaced)  =  43 

echo , 

echo,  Execute  after  OCMPILOT 


echo,  Hit  function  key  FI  to  abort,  any  other  key  to  continue 

pause 

state 

echo, - compute  yp  ,  yc,  and  yp*yc 

frequency ,p24 ,41 ,42,43 
kill,$$yp  4  rename, $$p24.g,$$yp 
p500=p(l ,1) 

frequency ,p500 ,41 ,42,43 
kill,$$yc  4  rename , $$p500 . g , $$yc 
data 

ypyc=yp*yc 

echo, - compute  phi 

state 

p500=p25( ,1:3) 

p501=diag(pl3) 

frequency ,p500 ,41,42,43 

kill,$$phi  4  rename ,$$p500 ,g,$$phi 

data 

phi=phi*p501*phi ’ 

phi=phi/p20(l , 1) 

phi=phi~ . 5 

phi=real(phi) 

kill , $$p500  4  kill , $$p501 

echo, - plot  result 

plot ,yp ,lwdm, all , auto, -60 ,60,6, "Y_p ,  Y_c,  Y_p*Y_c,  and  Phi" 

A 

YC 

ALL 

2 

A 

YPYC 

ALL 

3 
A 

PHI 

ALL 

4 
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OCMFREQ2 


echo,  0CMFREQ2  Macro 

echo,  Purpose:  Create  freq  plot  of  Y_p,  Y_c,  Y_p*Y_c,  and  Phi 
echo,  Uses  exact  calculation  of  delay 
echo , 

echo,  tau  =  Stl 
echo,  low  freq  =  &2 
echo,  high  freq  =  &3 
echo,  #  points  =  St4 
echo , 

echo,  Execute  after  OCMLQR,  OCMSETUP,  and  DCMKBF 

echo,  Hit  function  key  FI  to  abort,  any  other  key  to  continue 

pause 

state 

echo, - Compute  Y_c 

p500=p ( 1 , 1) 

frequency ,p500 ,&2 ,&3 ,&4 
kill,$$yc  k  rename, $$p500.g,$$yc 

echo, - Compute  Y_p 

p5C0=(p8(a)-pl5*p8(c) ; -p9*pl5*p8(c) ,p8(a)-p8(b)*p6) 

p501=(pl5  ,p8(b) ;p9*pl5) 

p501=(p501 ,p501 ( , l)+p500*p501 ( ,2) ) 

p502= (0*p6 ,  -p6) 

p503= (0 , 0 ,0 ,p502*p501 ( ,2) / 

pack ,p500 ,p501 ,p502 ,p503 ,p500 

frequency ,p500 ,&2 ,&3 ,&4 

kill,$$t  k  rename , $$p500 . g , $$t 

int,p500  k  dig , a  1  &  f requency ,p500 ,42 ,&3 , &4  k  analog 

kill,$$delay  k  rename , $$p500 . g ,$$delay 

frequency ,p7 ,&2 ,&3 ,44 

kill,$$nm  k  rename, $$p7.g,$$nm 

p501=diag(pl3) 

'  t  =  (H , F , H_3 ) ,  delay  =  exp(-s*tau) ,  nm  =  neuro-muscular 
data 

yp=nm*delay I -t( ,3)*t( ,4) 
ypyc=yp*yc 

echo, - Compute  Phi 


delay=delay I -t ( ,3) 
delay=delay*t( ,1:2) 

delay=(delay , 1) 

) 

phi=yc | -yp 
phi=phi*nm 

phi=phi*delay 

} 

phi=phi*p501*phi ’ 
phi=phi/p20(l ,  1) 
phi=phi* . 5 
phi=real(phi) 

kill , $$p500  &  kill , $$p50l  &  kill,$$p502  &  kill,$$p503 

kill , $$delay  &  kill,$$nm 

echo, - plot  result 

plot.yp.lwdm.all, auto , -60 , 60 , 6 , "Y_p ,  Y.c,  Y_p*Y_c,  and  Phi" 
A 

YC 

ALL 

2 

A 

YPYC 

ALL 

3 

n 

PHI 

ALL 

4 
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echo,  Enter  controlled  system  (2  different  ways) 
state 

p40=(-2, 0,0,1;  1,0, 1,0;  0,1, 0,0;  1,0, 1,0) 

cc 

yc=l/s 

yw=l/(s+2) 


M 


®ocmyc2 ,yc ,yw,p40 
state 

p40=p40(s , (2 , 1) ) 
p40 

echo,  Solve  the  OCM 
®ocmlqr,p40,l, .00017, .08, .001 

Qocmsetup ,8 . 8 , .15,-20,-20,-25, .00371, .09687, .04315,0,0,1 
Socmkbf , . 1 
Qocmpilot , . 15,2 
state 

gep,p24,yp 

yp=-yp 

echo,  Lou  order  tf  approximations 
near ,yp  ,ypl , 1 , le-4 
If approx ,ypl  ,yp2 ,5 , 1 

EX2 

echo,  Enter  controlled  system 

yc=l/s*2 

butter , yu  ,  .5,2 

Socmycl ,yc  ,yw,p40 

echo,  Solve  the  OCM 

Qocmlqr ,p40 , 1 ,6 .4053e-5 , . 1 , . 001 

Qocmsetup, 1 , .2 , -20 . -20 , -20 ,5 . 1223e-4 ,3 . 7070e-3 ,8 . 3460e-2 ,0 ,0 , 1 
Oocmkbf , . 1 

Echo,  Low  order  tf  approximations 
state 

gep ,p24 , yp 

cc 

yp=-yp 

near ,yp,ypl ,  1 , le-4 
If a , ypl , yp2 ,5,1 
near , yp2 , yp3 , 1 , .2 


S2 


EX3 


« 
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echo,  Enter  controlled  system 

senter ,yc ,3 ,0,4 , 1, .04,1, .9, 3, 1,0, 2, .  7 , .25,1,5 

senter ,yw, 1 ,0, . 2219 , 1 ,2 , . 7 , . 5 

Socmycl ,yc ,yw,p40 

echo,  Solve  the  OCM,  full  attention  case 
Socmlqr ,p40, 1 , . 00018 , . 1 , .001 

Qocmsetup , 1 , .2, -20, -20, -20, .0002707, .002026, .006276, .015, .025,1 
Oocmkbf , . 1 

Echo,  Low  order  tf  approximations 

Qocmpilot ,.2,2 

state 

gep ,p24,yp 

cc 


yp=-yp 

near ,yp ,ypl , 1 , le-4 
lfa.ypl ,yp2 ,8,1 
near ,yp2 , yp3 , 1 , .3 
Echo,  Fractional  attention  cases 
p 1 2 ( 1 , 6) =1/2  a  Qocmkbf , .1 
pl2(l,6)=l/3  &  Qocmkbf , . 1 
pl2(l,6)=l/4  a  Oocmkbf , . 1 
pl2(l,6)  =  l/5  a  Qocmkbf  ,  . 1 
pl2(l,6)=l/6  a  Qocmkbf,.! 
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